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Foreword 

Compiled in this volume are seven papers from agencies working 
with the Guidance Laboratory of NASA/ERC. These papers are of special 
studies in the disciplines of trajectory analysis, astrodynamics , and 
celestial mechanics. 

They include: 

(1) An extension of variational theory to cover problems 
involving functions that can be represented approximately 
only, through as closely as desired; 

(2) A presentation of an orthonormalization procedure for 
achieving a least squares approximation of a multivariate 
function; 

(3) A development of an analytic solution for minimum fuel 
impulsive transfer between low eccentricity orbits; 

(4) A development of a method for calculating ’’geometric" 
bounds for conditions where, if the total energy of the 
three-body problem is negative, then one body will recede 
to infinity from the other two bodies; 

(5) A development of a power series for the problem of three 
bodies where the coefficients in the series are generated 
by reversive operations; 

(6) A development of rigorous error bounds for approximate 
satellite orbit theories; 

(7) A development of a general perturbation theory for the 
long-term behavior of high eccentricity orbits about Mars. 

The first paper, along with extensions to this work, will be useful in 
trajectory analysis and guidance theory. The second paper should support 
trajectory analysis and guidance theory in supplying approximations to 
functions where only numerical values are available. The third paper will 
contribute to mission design and, in general, to astrodynamics studies. 

The fourth paper will be useful for trajectory analysis on mission design, 
along with contributing to studies in celestial mechanics. The last three 
papers will be useful in celestial mechanics, specifically orbit determi- 
nation or prediction. 



Guidance Laboratory Guidance Theory and 

Trajectory Analysis Branch 
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SUMMARY 


This volume contains technical papers on NASA- sponsored 
studies in the areas of trajectory analysis and guidance 
theory. The studies are being carried on by several uni- 
versities and industrial companies. These papers cover a 
period ending October 1, 1966. The technical supervision of 
the contracts is under the personnel of the Guidance Labora 
tory . 



Contents 


INTRODUCTION 

W. E. Miner, GST, ERC 

FUNCTIONS OF RELAXED CONTROLS 

Dr. J. Warga, Northeastern University 

DEVELOPMENT OF MULTIVARIATE FUNCTIONAL MODELS 

BY LEAST SQUARES 

Dr. D. E. Dupree and Dr. R. L. Truax, 

Northeast Louisiana State College 

A GENERAL SOLUTION FOR MINIMUM IMPULSE TRANSFERS 

IN THE NEAR VICINITY OF A CIRCULAR ORBIT 

Mr. T. N. Edelbaum, Analytical Mechanics 
Associates, Inc. 

REJECTION TO INFINITY IN THE PROBLEM OF THREE 

BODIES WHEN THE TOTAL ENERGY IS NEGATIVE 

Dr. D. C. Lewis, Control Research Associates 

AN ALGORITHM TO OBTAIN SERIES EXPANSIONS IN THE 

THREE-BODY PROBLEM 

Dr. P. Sconzo, IBM 

RIGOROUS ERROR BOUNDS ON POSITION AND VELOCITY 

IN SATELLITE ORBIT THEORIES 

Dr. J. V. Breakwell and Dr. J. Vagners, 

Stanford University 

AN INVESTIGATION OF HIGH ECCENTRICITY ORBITS 

ABOUT MARS 

Dr. J. V. Breakwell and Mr. R. D. Hensley, 
Stanford University 


Page 

i ojr 

7 i/' 


21 i/ 


37 




73 l/ 


121 





v 



Introduction 


Compiled in this volume are seven technical papers from 
agencies working under contract, or grant, to NASA's Elec- 
tronics Research Center (ERC) in the fields of guidance 
theory and trajectory analysis. This work was sponsored by 
the Guidance Laboratory at the NASA Electronics Research 
Center . 

The following table presents the contributing institution, 
the section of greatest relationship, and the discipline of 
the paper . 


SECTION INSTITUTION/ COMPANY 


Tra j ec t ory 
Tra j ec t ory 
Astrodynamics 
Celestial Mechanics 
Celestial Mechanics 
Celestial Mechanics 
Celestial Mechanics 


Northeastern Univ. 

Northeast La. State 

AMA 

CRA 

IBM 

Stanford Univ. 
Stanford Univ. 


DISCIPLINE 

Calculus of Variation 
Functional Models 
Impulse Transfers 
Celestial Mechanics 
Celestial Mechanics 
Celestial Mechanics 
Celestial Mechanics 


The present organization for the 
shown in Figure 1. 


laboratory effort is 


The following are reviews of 


the individual papers. 


Paper No. 1 

The first paper by J. Warga of Northeastern University 
is a contribution toward the development of variational theory 
by the methods of modern analysis. Applications are expected to 
be made on practical problems that cannot be handled satisfac- 
torily by the standard variational theory. The basic technique 
is to imbed the set of admissible control functions in a larger 
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•COLLFGE WORK STUDY PROGRAM 


FIGURE 1 

set that possesses an extremizing solution for the problem, 
and then to approximate that solution as closely as is desired 
with one of the original control functions. The present paper 
presents and proves the required existence and approximation 
theorems, and extends previous work on this theory done by 
Prof. Warga. Future work will include the development of the 
corresponding necessary conditions for an extremum and the 
consideration of constructive methods by which the solutions 
may be determined in practice. 

Paper No. 2 

The second paper is a technical progress report from 
Northeast La. State by D. E. Dupree and R. T. Truax. This 
paper first reviews the development of an orthonormalization 
procedure for achieving a least squares approximation of a 
multivariate function. The authors then go through the develop- 
ment of a weighting function to augment the least squares 
approximation. This is done in an attempt to require the 
resultant error to be less than a specified error tolerance. 

As yet the work is incomplete in that the proof only succeeds 
in showing that the weighted error will not exceed that of the 
unweighted function. 

Multivariate functional models are useful whenever func- 
tional approximation is needed for numerical data. Examples 
are cut-off requirements, time-of-ignition , and steering angles 
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which can often be numerically, but not analytically, generated. 
The objective here was to control the errors of such approxi- 
mations. Work will continue in this area. The basic ortho- 
normalization procedure has been programmed at ERC and will be 
used in application. To date no work has been done in imple- 
menting the weighting function. 

Paper No. 3 


The third technical paper, prepared by T. N. Edelbaum of 
AMA, entitled "A General Solution for Minimum Impulse Transfers 
in the Near Vicinity of a Circular Orbit," develops an analytic 
solution for the minimum fuel transfer between low eccentricity 
orbits. The author shows that in all cases a two-impulse 
transfer suffices. The results are derived by applying Lawden's 
primer concept. The range of validity of the results is deter- 
mined by the applicability of the linearized equations of 
motion. 

The fact that analytic solutions of the two-point boundary 
value problem were obtained makes this paper a valuable con- 
tribution to the state of the art in trajectory optimization 
and guidance theory. The impulsive solutions are ready approx- 
imations for finite time trajectory solutions and may be used 
as such in guidance modes for on-board as well as ground-based 
applications. 

Paper No . 4 

The fourth technical paper, on Rejection to Infinity by 
D. C. Lewis of Control Research Associates, sharpens a result 

the three— body problem due to Birkhoff . Roughly speaking, 
the result is that if, for given total energy (negative) and 
total angular momentum, the three bodies are initially suf- 
ficiently close together, then one of the three will ulti- 
mately recede infinitely far from the other two. The result 
is, however, a very precise one, and a rigorous proof is 
given. The condition on the initial configuration is stated 
in terms of the Lagrange inertial radius R which is a measure 
of "closeness together" of the three bodies. A method is 
developed for calculating three quantities R q , and t ± such 

that if initially R >_ R^, then r at t^, R > R^. Further, the dis- 
tance between the two bodies closest together at will remain 
bounded (with an estimate given for the bound) thereafter, and 
the third body will recede to infinity from these two. 

This study was done to get a better understanding of the 
problem of three bodies in celestial mechanics. Establishing 
the bounds is an advancement which allows us in-house to attempt 
to develop methods of studying analytically swing-by trajecto- 
ries to distant places. 
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Paper No. 5 

The fifth technical paper, entitled "An Algorithm to 
Obtain Series Expansions for the Three-Body Problem," by 
p Sconzo, IBM, deals with power series solution, both in 
the time variable itself and also in terms of a regularizing 
variable. Outlines are presented for the recursive generation 
of the coefficients in the two series. A lower bound is given 
for the radius of convergence of the series in time and the 
relation of this series to the "f" and "g" series representa- 
tion of the solution of the two-body problem is briefly ais 
cussed. (It should be noted that the method is of wider 
applicability and, in particular, the relation with the f 
and "g" series might well be exploited for artificial satellite 
theory.) Next the properties of a general class of regulariz 
ing transformations are presented with special reference to t e 
transformation of Levi-Civita. Subsequent reports will deal 
with convergence of the series in the regularizing variable 
and application of a machine-generated series based on the 
recursive formulas mentioned above to a problem of Bohlen for 
which numerical calculations by Zumkley are available for 
comparison. Also included are some comments on a solution of 
the three-body problem recently obtained by Bazayevesky in 
terms of a power series in time. 

While the method used in this report has been previously 
applied (Steffenson, Fehlberg, Rabe et al) , there are some 
novel features in its development and in the symbolic program 
for the generation of the series. The theory developed in 
this report is directly applicable to problems for which the 
three-body problem forms a good model, for example, to ephemens 
calculations for motion about the sun of two planets or one 
planet and its satellite. 

Paper No. 6 

The sixth technical paper by J. V. Breakwell and 
j. Vagners discusses rigorous error bounds for approximate 
satellite orbit theories. The paper is concerned with the 
error in prediction from initial conditions over a time inter 
val of order 1/e for a general perturbation theory developed 
in powers of e. The problem is that if one, or more, of the 
variables possesses a linear dependence on time, the coeffi 
cients of time in these variables must be calculated to second 
order in e in order to obtain a first-order approximation for 
time intervals of order 1/e. The development is carried out 
for earth satellites (e = J 2 ) using the Brouwer-von Zeipel 
technique with Poincard variables. The only one of these van 
ables with a linear second-order (in e) time dependence is an 
anqle closely related to the mean anomaly. The authors make 
use of the energy integral to obtain this second-order coeffi- 
cient, thus bypassing the necessity of a full second order 
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formulation and integration of the perturbation equation. The 
theory includes tesseral harmonics, but omits consideration of 
resonance effects arising from the tesseral harmonics and also 
from critical angles of inclination. These will be the object 
of future investigations. The theory is applied to a circular 
2000-mi satellite for 7 days, and the resulting errors are 
consistent with the theory. 

Paper No, 7 

The seventh technical paper, by J. V. Breakwell and R. 

D. Hensley of Stanford University, develops a general per- 
turbation theory for the long-term behavior of high eccen- 
tricity orbits about Mars. Mission requirements on a planetary 
orbiter would probably require a small pericenter distance 
(for observational purposes) and high eccentricity for saving 
of fuel. In this study, a pericenter requirement of 4000 to 
6000 km was imposed and eccentricity larger than 0.5 was taken. 
For such an orbit about Mars, the orbital period is small com- 
pared to the Martian year which, in turn, is small compared 
to the rates of change in the orbital parameters due to Mars' 
oblateness and the sun. These facts make a double averaging 
procedure, first, over the orbital period and then over the 
Martian year useful for the study of long period effects. The 
dominant perturbation is due to Mars' polar oblateness and 
results in secular rates for the argument of the node. The 
perturbations caused by the sun result primarily in long period 
fluctuations in inclination and eccentricity and hence in peri- 
center distance. In addition to the oblateness critical angle 
of 63.4 degrees, a number of other critical angles occurs. It 
is in the neighborhood of these critical angles that the ampli- 
tude of the fluctuations in eccentricity are most pronounced. 

A detailed analysis of these resonance effects is given. In 
addition, an appendix contains an analysis for small eccentric- 
ity (J^ for Mars). As a first step in the analysis of plane- 
tary orbiters, this is a very useful study and extension to 
include short period effects would be desirable. The applic- 
ability of the analysis to orbiters of other planets would 
require consideration of the relative magnitudes of the orbital 
period, the planet's period about the sun, and the rates of 
change of the orbital elements, as well as a study of which 
perturbations are dominant. Development of theories for 
planetary orbiters is, of course, essential to any extensive 
program of planetary exploration. 
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Functions of Relaxed Controls* 


by J. Warga 
Northeastern University 
Boston , Massachusetts 


summary N67-29371 

Problems of the calculus of variations do not admit, in general, ordinary 
minimizing solutions. General existence theorems in the calculus of varia- 
tions have been established only with the introduction of generalized, or 
relaxed, solutions of the given differential equations. These relaxed solu- 
tions represent the limits of ordinary solutions with rapidly oscillating 
derivatives. 

In the present paper we consider a class of problems of the mathematical 
control theory that are not necessarily defined by systems of ordinary dif- 
ferential equations; they may involve solutions of certain partial differential 
equations, nonadditive set functions, or other functionals. The controls are 
functions p from a metric space T, with a given measure, to a metric space 
and are subject to restrictions of the form p(t)eR#(t) a. e. in T, where 
Rff (t) is, for every t in T, a given nonempty subset of T. We consider 
functionals x(p, b) = (x*(p, b), . x n (p, b) ) depending on controls p and on 
parameters b restricted to a compact set Bq. The variational problem con- 
sists in determining the minimum of x^(p,b), subject to the previously 
mentioned restrictions on p and b and the condition that x(p, b)eBj, where 
Bj is a given set in the euclidean n- space. 

We establish an existence theorem asserting that, in a large class of 
such problems, the restricted minimum is achieved by n relaxed" controls. 

We also establish an approximation theorem stating that each such relaxed 
control can be constructively approximated by original, or ordinary, 
controls. 

Necessary conditions for minimum will be discussed in a forthcoming 
report. 

* This research was initiated at Avco Corporation, Wilmington, Massa- 
chusetts, under N. A. S. A. -ERC contract ^AS—12- 1 12, and has been con- 
tinued at Northeastern University under gfrant NGR 22-011-020, 
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1. Introduction. The mathematical control theory deals primarily with 
functionals defined in terms of a system of ordinary differential equations. 

It is our purpose here, and in a paper to follow, to extend certain methods 
and results of the control theory to a more general setting. In particular, 
we wish to generalize certain results of references [1] and [2], In this en- 
deavor, and especially in arguments pertaining to problems of existence, we 
continue to apply certain concepts first introduced by L. C. Young in his 
study of M generalized curves" [3], [4], 

Let 5? be a class of mappings from a set T to a set R, and let x (i=l, ,.,n) 
be real-valued functionals on ft , that is, functions from 51 to the real line. 

In many problems, the sets T and R are metric and the vector x of function- 
als is characterized by a "weak" continuity, in the sense that x(p) and x(p) 
differ little if pe ft , p"e ft and p (t) and p (t) are at a small distance from one 
another for all (or almost all) values of t in T. This type of continuity (with 
respect to uniform convergence) has been most frequently investigated in 
connection with differential equations (e.g. in the study of perturbation 
methods), in the calculus of variations (weak variations), etc. 

In a large class of problems, the functionals x are also continuous in a 
different sense. Let us consider, as an example, the problem of determining 
the temperature 0 (t^) at a time t and at a fixed point ^ in the interior of 
a conducting body whose surface is subjected to a heat flux that varies with 
time and position. If the heat flux is only slightly changed at all times and 
over the entire surface, the value 0 (t, £) will change but little. This is 
so because 6 (t, f ) is a continuous functional of the heat flux in the pre- 
viously described sense. Assume now that the interval of time [ 0, t ] is 
subdivided into 2k equal subintervals, and that at every point z of the surface, 
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the flux h(t, z) during the interval [ t, ~ t” ] is replaced by 
h(t + , z) and the flux h(t, z) during the interval [ T, ? t~ ] 

is replaced by h(t- z). for i = 0, 1, . . , k- 1. Then for large values of k, 
we would expect 0{t, to be affected little by this "permutation", even if 
the flux h(t, z) is a rapidly varying, or even a discontinuous, function of t. 

We would also expect 6 (t, £) to be little affected by " permuting' ' the heat flux 
between many adjacent small areas of the surface of comparable measures. 
This second type of continuity is relative to a mode of convergence resem- 
bling the weak convergence of measures. It is of fundamental importance in 
control theory and in our further considerations. 

The continuity of functionals in this second sense makes it possible to 
simulate mathematically the limits of rapidly oscillating functions in 91 . 

We shall refer to elements of !R (functions from T to R) as "original con- 
trols", and we shall imbed $ in a larger space S of " relaxed controls." If 
we assume that R is a compact Hausdorff space, then we can define the 
class S of regular probability measures on Borel subsets of R. A ” relaxed 
control" a is a function from T to S. The relaxed control a will simulate the 
limit of rapidly oscillating ordinary controls p j, p^, .. when for all (or 
almost all) t and all Borel subsets R^ of R, the a (t) -measure of Rj re- 
presents, in some sense, the limit, as j -*oo, of the relative frequencies of 
occurrence inside Rj of the points pj(r) in the neighborhood of t. A relaxed 
control <j with the property that the measure o(t) is, for every t, concentra- 
ted at a single point p (t), can be identified with the original control p . 

Relaxed controls, patterned after L. C. Young' s definition of "general- 
ized curves", provide a means of completing the space $ of original controls. 
Their use paves the way, in complete analogy with the calculus of variations, 
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to existence and approximation theorems which we state in section 2 and 
prove in section 3. The next task is the derivation of necessary conditions 
for minimum. We propose to discuss this subject in a paper to follow. 

2. Existence and approximation theorems. We shall assume in this 
section, and in section 3, that T and R are compact metric spaces, that 
B Q is a compact Hausdorff space, and Bj is a closed set in the euclidean 
n-space E n - We shall assume, further, that a nonnegative, finite, regular, 
complete, and nonatomic measure is defined on T. We represent by p, and 
sometimes by p ( ), a mapping from T to R, and by p (t) the image of a point 
t under the mapping. A similar distinction is consistently made between a 
function (mapping) and the image of a particular point under the mapping. A 
mapping from T to R is " continuous" at t (on a set Tj) if | p (t), P (f ) | — 0 
as t . t for all teTj). Here | r y r 2 | designates the distance of 

points rj and r 2 in R, and similarly | t' , t | will designate the distance m 
T. A mapping from T to R is "measurable" if, for every e > 0, there exists 
a closed set F £ in T, whose measure | F £ | is at least | T| -t, and such 
that the function p is continuous on F £ when restricted to F £ . 

These definitions of continuity and measurability of p can be easily 
seen to be equivalent to the following statement: p is continuous at t (on a 

set Tj), respectively measurable on Tj, if the function ip , defined by 
0(t) = >Mp(t)) on T, is continuous at t(on Tj), respectively measurable on 
T , for every choice of a continuous function <j) from R to Ej, 

Let R# be a mapping from T to the class of nonempty subsets of R, and 
let 3 ! be the space of measurable functions p from T to R. We are given a 
function x: K x B fl — E n - We wish to investigate the original problem of de- 
termining the minimum of x'(p, b), subject to the restrictions that 
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p(t)eR#(t) a. e. in T and x(p, b)eB^. Alternately, we wish to consider 

" approximate minimizing solutions” to this problem. An "approximate 

00 

solution" x is a sequence {p . b. } such that the p. are measurable 

J J=1 J 

functions from T to R, b-eB p.(t)eR#(t) a. e. in T, x(p . b.) converges, 

J u J J> J 

as j— oo, to a point x^, and x^cB^. An "approximate minimizing solution" 
is an approximate solution that minimizes . 

Let now S be the class of regular probability measures defined on the 
Borel subsets of R. We shall refer to a function p from T to R as an 
"original control", and to a function o from T to S as a " relaxed control". 

A relaxed control cr is "continuous", respectively " measurable" , on a set Tj 
if {dr;t) is a continuous, respectively measurable, function on Tj 

for every choice of a continuous function (f> : R-*-Ej. Here a(R^;t) repre- 
sents the cr{t) -measure of a Borel set Rj. We can easily verify that if a is 
a measurable control and Rj is a Borel subset of R, then<j(Rj;t) is 
measurable. We shall denote by S the set of measurable relaxed controls. 

If a relaxed control has the property that (t) is a measure con- 
sisting of a single mass point p (t) a. e. in T, then we refer to it, somewhat 
loosely but without any fear of confusion, as the original control p . In this 
sense we consider 91 to be a subset of S . We shall also treat original, 
respectively relaxed, controls as identical if they differ only on a set of 
measure 0 in T. 

Definition Z, 1. The Young representation. We shall say that a function 
y: S xBq— is a" Young representation of x" if y coincides with x 
on $ x Bq ; that is, if y( a , b) = x(p, b) for every b in Bq and every relaxed 
control cr ^ such that (t) is concentrated at the single point p (t) a. e. in T. 
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(We observe that the relaxed control is measurable if, and only if, 

the original control p is measurable). 

Let now L^T) be the Banach space of real -valued integrable functions 
on T, and let C(R) be the Banach space of real -valued continuous functions 

on R, both with their conventional norms. ( I C(R) = max l I ^ II (T) ~ 

/ , v 1 reR 1 

T | f | ) . We shall define functionals k(4>, f, a) on QRJxL^T) x § by 


k{0, f, a) 



(r)cr (dr;t). 


If a is an original control p, we shall write 
P 


k(</>, f, o p ) = k(<£>, 


» f » P ) *“ ^ n-> 


f(t)<Mp(t». 


Definition 2. 2. The Young topology on S . We shall say that a sequence 

CTil a of . . . in S is convergent if the sequence of real numbers {k(4>, f, a) } "j 
12 J J 

is convergent for every choice of ( <f> , f) in C(R)xL.^(T). We shall say that 
a in S is a limit of the sequence Oy Oy ... if 


k(i£, f, a) = lim k(4>, f, a.) on C(R) x L^(T). 
j-00 J 

# 

We now consider a mapping R satisfying 
Assumption 2.3, For every e >0 there exists a closed subset T £ of T, of 
measure at least | T | - e , with the property that 

(2.3. 1) for every t"e T £ and every r€^(F)(the closure of R#(t)) there 
exists a measurable original control p, continuous at t when restricted 

u 

to T £ , and such that |p(t)»r| <e and p(t)eR (t) on T; 

(2.3.2) the mapping R# , when restricted to T £ , is continuous with 
respect to inclusion, i. e. for every t in T £ and every h> 0, there 
exists 6 = 6(h, F) such that R#(t) C U(R # (t).h) and R # ( T) C U(R#(t), h) 
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if t e T £ and | t, t | < 6. Here 

U(R h) = {r e R| | r, r j | <h for some r j in Rj }. 

We can now state our basic approximation and existence theorems. The 
symbol R# (t) will denote the closure of R^(t). 

Theorem 2. 4 . Let the mapping R^ satisfy condition (2. 3. 1). Then every 

measurable relaxed control o can be approximated (in the Young topology) 

by measurable original controls P P £* • • • • ^ o(R^ (t);t) =1 a. e. in T 

u 

then the controls p . . . can be chosen so that p^(t)eR (t) a. e. in T 

(j = 1,2...). 

Theorem 2. 5. Let S ^ = {oe S | a(R#(t);t) =1 a. e. in T}, and assume that 
the mapping R^ satisfies Assumption 2. 3. Then the set § ^ is sequentially 
compact. 

Let y be a Young representation of x, and assume that y is con- 
tinuous on S ^xB q (with respect to the product topology on S ^xB Q ). Let 
ft # = {p€ ft | p(t) e R^(t) a.e. in T), X = {x(p, b) |pe ft be B Q } and 

Y = {y(a,b)| ae S beB 0 }. Then Y is the closure of X. 

As a corollary of Theorems 2.4 and 2. 5, we derive 

Theorem 2, 6, Let ft ^ and X be defined as in Theorem 2. 5, and let the 
assumptions be the same as in Theorem 2.5. Then either YPlBj is 
empty, or there exist a € S ^ and Td e Bq that yield the minimum of 
y^a, b) subject to the condition y(cr, b)eBj. If the construction described 
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^ 

in 3. 3 is used to approximate ct with a sequence p ^ p • in K , then 

the sequence {p\, b is a minimizing approximate solution of the 

original problem. 


3. Pr oofs of the approximation and existence theorems. 

Definition 3. 1. A dense sequence of partitions of T [5, pp. 171-174], 

We shall say that P T is a dense sequence of partitions of T if 
P T = {P T 1 ,P T 2 ,..}; = {T j. T^, • • > } (i = 1 , 2, . . the sets 

Tj (j = 1, . . , j.) are, for each i = 1, 2, . . , measurable and disjoint and 
jjji rpi _ every element of P T i+1 is contained in some element of 

j=l j T 

p^ f for i = 1, 2, . . ; and to every measurable subset E of T and every 
e>0 there correspond a positive integer i(e ) and a subset J(E, e) of 


i( £ ) 


{ 1.2 j i(c) } suchthat |E-E 0 | + |E 0 -E|<e. where E Q = U je J(E> £ ) T j • 

It is well known [ 5, Th. C, p. 173] that there exists a dense sequence of 
partitions of T as a consequence of T being metric and compact, and the 
measure on T having the properties listed at the beginning of section 2. 

We shall require a lemma. 

Lemma 3.2. Let e >0, T £ have the properties described in Assumption 
(2.3.1), F be a measurable subset of T £ , {R - . . $ a partition 

of R into disjoint nonempty Bor el subsets, and ae S # and assume that the 
support of a (t) is contained in R^ (t)(the closure of R#(t)) for all t in F. 

Let a k =f F a(R k ;t) (k = 1, . m). Then there exist a partition of F into dis- 

joint measurable sets F|, F^, . F^ and a measurable original control p 
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such that, for k = 1, 2, . . . , m, | F k | = a k , p (t) e R#(t) on T, and p(t) is within 
a distance 2e of a. e. in F^. 

» 

* Proof. Let k represent integers from 1 to m, and let 

G^= {teFj a(R k ;t)/0 }. For every nonempty subset A of { 1, 2, . . , m } , 

m 

let G^= Since ^ a (R^jt) = a (R; t) = 1 in T, it follows that 

k= 1 

F = G^ (the union over all nonempty subsets A of { 1, 2, . . , m } ). For 
every set A, we can partition G A into disjoint measurable subsets 

k f k 

G^(keA) of measure J q o-(R k ;t). If k?A we define G^ to be the 

empty set. We now let F^ = G^ , and verify that | F^J = a (k=l, . m). 

Let now k be fixed. Since, by construction, cr(R k ;t)^0 for teF k , for 

every t in F^ there exists a point r k in R^ (t)O^ and, by Assumption 

t k _ 

(2, 3. 1), there exists a measurable original control p_ , continuous at t 

i k k - * k u 

when restricted to F, and such that | p_(t), r_ | <e and p_ (t)eR ff (t) on T. 

k _ ^ 

Because p_ is continuous at t when restricted to F, there exists a 

* _ _ k k . 

neighborhood (relative to F^) N k (t) of t such that | p_(t), r_ |<2e in 

_ k _ t t 

N^ft); hence p_(t) is within 2e of R^ for tcN k {t). Since is covered by 

open neighborhoods (relative to F^) N k (t), it must be covered a. e. by a de- 

_ k 

numerable subfamily, say N k (tj), N^ft^), ... . We now let p (t) = p_ (t) 

_ 1 ^ 

fortcN,(t.) (k= 1, . . . , m; j= 1, 2, . . . ) and p (t) = p (t) everywhere else on F. 

k J F i 

Since p_ (t) e R^(t) for all k and j, it follows that p(t)eR^(t) on T. We also 

t j 

observe that p is measurable and p (t) is within a distance 2e from R^ a. e. 
in F k (k= 1, 2, . ., m). 

3. 3. Proof of Theorem 2. 4. Let the sets (j=l, . . . , j.;i=l, 2, . . ) define a 
dense sequence of partitions of T as in Definition 3.1. Since R is metric 


15 



TRAJECTORY ANALYSIS AND GUIDANCE THEORY 


and compact, for every positive integer i we can partition R into disjoint 

Borel subsets R^ (k= 1, . . , k^ ;i= 1, 2, . . ) of diameters not exceeding i . In 

every one of these sets R^ we may arbitrarily select a point rj^ . Let 

Tj/i be defined as in Assumption 2. 3, and let T* 1 = T^ OTj^. 

(j=l, 2, . . , j.;i=l, 2, . . ). For every fixed positive integer i and for every 

j= 1, 2, . . . , j., we may define sets T* , (k= 1, . k.) and original controls 

l j, k l 

p* that have the properties described in the statement of Lemma 3.2, with 

/ #i i i i 

1/i, T. , T. ,, R, , and p. replacing £, F, F, , R, , and p , respectively. 

J J» K K J K K 

Let now a measurable original control p^ be defined for i=l, 2, . . . by the 

relations . 

P ^t) = p j{t) on T. 1 (j= 1,2,.., j.), 

P^t) = P \(t) on T-T^. 

We observe that p . (t)eR^(t) on T, p . (t) is within a distance 3/i of r* a. e. 

in k . and | T‘ ( k | = J T * j ct (Rj.;t) (k= 1. 2, . . , k.) . 

Let now e > 0, <£cC(R), and let E be a measurable subset of T. The 

symbol I (j> | will denote the C(R)-norm of <f>, i. e., Max \ <f){r) | . We may 

re R 

choose an integer i^ sufficiently large so that, for every i^i^, there exist 

a subset J. of { 1, 2, . . . , j. } and a measurable set E. in T such that 
i 1 J l 

E i = jVj. T j' 1 * | E-E.| + | E.-e|< I E /|*|and|*(r)-*<r')|<:J c/|t| if|r, r'|<3/i. 

Finally, let 0(a) represent here a quantity not exceeding a in absolute 
value. Then, for all i> i^, 
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= E 

jej. 


k. 

1 


E 


k= 1 





k i f 

= E E J j #(p i (t» + o(E) 

jeJ. k=l T* 1 4 

1 J, k 


*,f, *<»!«>+ »<!■> 

/ J 

= •/ £ * (PjW) + 0(c). 


Thus 



^eC(R) and every measurable characteristic function f_. It follows that 

E 

lirn^k^, f,p .) = k(*. f, CT) for every (0, f) eCfRJxL^T). This completes the 
proof of the theorem, 

The proof of Theorem 2. 5 is largely a generalization of a construction 
of L. C. Young [3], 


3 ' 4 El°°l of Theorem 2. 5. We shall first prove that S is sequentially com- 
pact in the Young topology. Let a j, cr^ . . . be a sequence in S , and let 

kj(^f) = k(*. f, a.) = J T f(t) J R 0( r)o.(dr;t) (*eC(R), feLj(T), J - 1, 2, . . .). 

The bilinear functionals k^ clearly satisfy the relation 

(3.4.1) |k.(<M)|<|f| |*| (*e CtRKfeLjfT), j=l,2,...), where 

l f l = I £ Il,(T) = f T | f W| an( M</> | = |* | = Max |*(r)|. 

1 ' re R 

Let now C* (R) be a dense denumerable subset of the separable space 

C(R). For all *eC' (R) such that | *|<1, the sequence of linear functionals 
00 

{kj($, )} on Lj(T) has norms bounded by 1, as a consequence of relation 
(3.4. 1). It follows that there exists a subsequence of the functionals k.(<£, ), 
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which we shall continue to designate by {k^(0, )}j = j, that converges to a 

bounded linear functional k(0, ) for all 0eC' (R). Relation (3.4. 1) clearly 

continues to apply to k(0, f) for all 0eC' (R) and feL^T). We can easily ex- 
tend, by continuity, the definition of k(0, ), as a bounded linear functional 

on Lj(T), for all 0eC(R) and we verify that k is bilinear and satisfies the 
inequality (3. 4. 1). 

For every fixed <f> in C(R), k (<f>, ) is a bounded linear functional on 
LjfT) and, as such, can be represented by 


L 


K t)f(t). 


where £(</>, t) is a bounded measurable function on T. Since, for each 
0, ?{0, t) can be arbitrarily changed on a subset of T of measure 0, we can 
determine a subset T 1 of T such that | T 1 j = j T | and £ ( , t) is a bounded , 
linear functional on C' (R) for all teT 1 . We then verify that F { , t) can be ex- 
tended, by continuity, for each teT*, to a bounded linear functional on C(R), 

and that /* 

k(0, f) = 7 T f (<*>. t)f(t) <0eC(R), feL^t)). 


Furthermore, relation (3.4. 1), applied to k(0, f), implies that there exists 
a subset T* of T' , of measure |t|, such that 
I t (+.t)|s|* | on C(R) xT*. 

We can prove this last relation first for all 0eC' (R), and then, by con- 
tinuity, for all 0 e C(R) . 

We can, therefore, conclude that there exists a signed regular measure 
cr(t) for all teT* such that 

(3.4.2) £ (<f), t) = J R 0(r)cr(dr;t) on C(R)xT*. 
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Thus 


lim k(</>, f, a.) = lim k.(<£, f) = k(4, f) = 
j— oo J j— oo J 



4(r)<j(dr;t) 


= k{4 # f,a) (feL^T), 4cC(R)). 


Since k.(4, f)>0 for all j if <j)( r)>0 on R and f(t)>0 on T, it follows that 
k(0, f) has this property and, therefore, £ (<f>, t)>:0 a. e. in T*, say in T*(<f>), 
if <J(r)>0 onR. Furthermore, k^(4j» f) ■Jr f(t) if 0j{r) = 1 on R, hence 
k (0j» f) = J ^ f(t), and it follows that i (0 Jf t) = 1 a. e. in T*. Since C(R) is 
separable, we can prove, as in previous arguments, that there exists a sub- 
set of T of T of measure j T J such that £ {<j> t t)2r0 on T^ for every non- 
negative 4 in C(R) ; and t) = 1 on T # if ^(r) = 1 on R. It follows then 

from (3. 4, 2) that a (t) is a regular probability measure for every t in T^. 
The corresponding mapping a is measurable on T since £ (<j), t) is measur- 
able for every <£eC(R). This shows that ae S and completes the proof that 
S is sequentially compact in the Young topology. 

We shall next show that if a sequence cr^, ct^, . . . converges to a in the 
Young topology, and if a(R# (t);t) = 1 a. e. in T (j=l, 2, . . . ) then cr(R#(t);t) = 

1 a. e. in T. Indeed, let rj>0 and let T^ of measure at least |T | - rj be a 
closed subset of T such that R^ is continuous when restricted to T . Let 
e>0, teT^, (t, 6) = {teT^ | | t, t |<6} f U(Rj, h) be the open h-neighborhood 

of a set RjCR, ^l/2 = U(R#(t), ^ 1/2 t ^ ie c ^ osure ^1/2* 

Uj = U(fT^ (t),e ), and let 6 = 6 (e) be such that R#(t) (T an< ^ 

R^( t ) C U(R^( t ), -=) for all t in S (t , 6) . Let in C(R) be such that 
<F(r) = 0 on tJ^ 2 , 0<^(r)<l on R, and ~${r) = 1 on R- XJy Then 
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0 = lim / / <?(r)<J (dr;t) = / __ J 

i . c /T" c \ D ^ f M R 


(r)tj(dr ;t) 


j— » S (T,6) 


S (t, 6 ) R 

n 


>/ a(R-U;t); 

J S ri (f.6) 


hence a(R- U 1 ;t) = 0a.e. in S (1,6(0) or a(Uj;t) = 1 a. e. in S^t, 6 (e)). 

We observe that, for all t in S^(t, 6(e)), U j C U(R ^(t), ^ )• It follows that 
T can be covered by open (relative to T ) neighborhoods in each of which 

n v 

a(t) is a. e. supported by U(R^(t), ). Since e is arbitrary, R (t) is compact 

jj. 

for all t, and the measure a (t) is regular, we conclude that a(R ff (t);t) = 1 a. e. 
in T . Since n is arbitrary, it follows that o(R^(t);t) = 1 a. e. in T. Thus 

n 

S ^ is sequentially compact. 

It follows from the above conclusion and from the continuity of y on the 
sequentially compact space S^x-Bq that the set Y is closed. By Theorem 
2. 4, Y is contained in the closure of X. Since X is obviously a subset of Y, 
we conclude that T is the closure of X. This completes the proof of the 


theorem. 
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Development of Multivariate Functional 
Models by Least Squares* 


by D. E. Dupree and R. L. Truax 
Northeast Louisiana State College 
Monroe, Louisiana 


N67-29372 

SUMMARY 4 " 7 ^ 

A technique for deriving an approximating function yielding an error, 
in the sense of least squares, less than a specified error tolerance is 
developed. 


A COMPUTATIONAL TECHNIQUE FOR DERIVING 
THE LEAST SQUARES APPROXIMATING FUNCTION 

Given : n + 1 tabular points {p o ,X(p Q )}, (^ 1 ,X(^ 1 )}, ... , (£ n ,X(£ n )} 

for the function X = X(f3), where £ = (x , x^, . . . , x^). 


Problem : Choose N + 1 independent functions q^(p), cp^(P )>•••; cp^(P) and 

determine the polynomial E A.<p.(p) satisfying the property 

j-0 J J 


that 

F(A 0 ,A i , 


>0 = E tx(p.) - E A.cp.(p.)} 2 
* i=0 1=0 J J 


is minimum. 

dF dF 

A necessary condition for F to be minimum is for = 

3A 0 SA 1 

i^F 

... = — — = 0. This yields the following system of 

3a n 

★Performed under NASA Grant NGR-19-006-001 
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normal equations : 

+ + ••• + V'Ph’V = »*> 

A^qtyq^) + A^q^q^) + ... + ^((fy = 


VV V + + ••• + " A^,) 

where X *= {X((3 q), X(p^)^ ... , X(p^)}> 

cpj * (cPj(^q)^ <pj(P^)j ••• > <pj(Pn^> ^ = 0,1,...,^ 

and <^) denotes the inner product E 

Using the vectors cp^j • • • > cp^j define a set of vectors 

VV* 4,, ® N as follows: 



This is the orthonormal collection yielded by the Gram- Schmidt 
orthogonalization process; that is ( e ^ e j) = ^ ^ j, and 

( e .*> e .s) = J = 0,1; 

o «J 

In addition, define the triangular array of coefficients 
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V- 1 ) 


A 1 (-D 

\(o) 




Ag( - l) 

Ag(0) 

V 3 -) 



Aj("l) 

*3(0) 

AjCl) 

*3(2) 


V- 1 ) 

v°> 

V 1 ) 

A n (2) -- 

A^N-l) 


where 

A (-1) = 1 

Y - - Y -1 _ _ 1/2 

{( V'Py ) - ^ ( V 6 0 ) } V 

A y (k) = , k= 0,1 ,..., y-1. 

«vV -j VvV ,1/a 

Then e, and A.(k) ; 0 < j < N and k = 0, 1, 2 , . . . , j^l, can be 
J 0 

written as follows: 

k-1 

Aj(k) = A^-^y-lKcpj,^) - E A 3 (i\W 

j-1 

e 4 = A (-l)tp - E A (i)e . 

J j j i=0 " ± 


Then the coefficients of the triangular array and the e^, 
0 < j < N, can be written recursively as follows: 
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To see how the coefficients of the triangular array and the 
"e , 0 < j < N, are used, define N + 1 functions as follows: 

J - “ 

f 0 o) = y-n^o) 

3 - 1 

Tj(p) = A^-lJcpjO) - S AjCiJf^p), 1 < J < N. 


Then each f (f3) is a linear combination of the N + 1 functions 
J 

cp 0 (f3)> cp^CP); ••• t and 

V {f J (f V' •- ’ V P n )) = V 

n N 

Thus, F(A q ,A 1 , {X(p i ) -2 AjtpjOj^)) 2 

n N 

= Z (X(p ) - Z A’f (P,)} s 

i=0 j=0 J J 


= F' (A^j, A^, • • • , Ajj) 


and the necessary condition that F' be minimum yields the normal 
equations 
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N N 

This yields the function E A'f.(p) = S A.m,(f5) such that 

j J=0 J J 

n N 

2 (X(p.) - 2 A.cp.(p. )} 3 is minimum. 
i=0 1 j=0 5 J 1 


DERIVATION OF A FUNCTION SATISFYING A GIVEN 
ERROR TOLERANCE IN THE SENSE OF LEAST SQUARES 


n N 

Although E = 2 (X( P ) - 2 A.<p (p )} 3 is minimized by the least 

i=0 j=0 J J 1 

squares procedure, there is no assurance of the relative size of this error. 

ir* 1 ?' We .^ e f d x ^° be able to determ ine an approximating function in such a 
fashion that the error, in the sense of least squares, will not exceed a 
given tolerance. Before doing this, notice that 


E 


n 

S 

i=0 




H x " A d'Po ‘ ” • • " A n'P N II 2 


= II X 


s (X,e )e || 3 
j=0 3 3 


II X II 2 


If 

[X, 2 (X, e )e ] 

>0 J J 


N 

- [X, 2 (X, e )e ] + 

J=0 J J 


II 2 (X,e )e I) 3 
J=0 J J 


- II X IP 


N 

2 (X,e ) 3 . 

J=0 J 


From this we note the following points : 

1« || X |p is an upper bound of E. 

2. A sum of any k of the N + 1 terms (X,!.) 3 , 0 < k < N, will yield an 
error E 1 > E. 

3* If e N+1 is any other non-zero vector orthonormal to each of 
- - - - N+ 1 

V e i'-**>V then II X II 3 - 2 (x,e ) a < E. 

j=0 0 
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Thus, 

Problem; 


Solution; 


we may define the problem stated above as follows: 

N 

After evaluating || X || 2 - £ (X, e ) 3 , we find that this 


value still exceeds a given error tolerance 6. Then we need 
to find such that 

II X II 3 - S (X,e ) 3 - (X,e N+1 ) a < 6, or 


i=o 


N 


(X, e N+1 )2 > || X || 2 - (X, ej ) 2 - 6, 


where e^ +1 is the vector associated with that is orthonormal 

to e 0 ^ * * * > e N « 


^ e t ^+1 “ (v x i^*-> x n ) # T ^ ien 


'PN+l " %+l " * ^ c PN+l ,e j^ e j 


J=0 


N n _ 

“ (^<y X J ' ' ' *' X n^ " ^ ( 2 

u x n , 0 i=Q i jj. j 


where = ( W * # e jn } ' ° - J - N * 


Then 


N+l 


N n 

u l n j =0 i=o J J 


n N n 


t s x? - s ( s V-ii) 3 ) 

i=0 1 j=0 i=0 1 Ji 


1/2 
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and if X = ( , t^), then 


N n 


[ 2 \.t - Z ( Z \ .e )(X,e )J 

* r\ Jw -i- i r\ i a 1 J 1 J 


( X ^N + l) 3 


i=0 


j =0 1=0 


N n 


S ^ - 2 ( S X^) 3 


i=0 j=0 i=0 


Thus, to have (X,e ) 3 > || X || 3 - E (X,e ) 3 - 6, 


j=0 


we must have 


n N n _ 

[ 2 U, - Z ( Z \,e )(X,e )] a > 

A — r\ x A a r\ a n 1 J- 1 - .J 


i=0 


j=0 i=0 


n N n 

12 


[ Z X* - 2 ( 2 X e ) 3 ] [ || X || 3 - s (X,e ) 3 - 6], or, 

1=0 1 0=0 i=0 1 j=0 0 


equivalently, we must have 


& *■ ( X,e 0^ e 0i "• ^ X,e N^ e Ni^ 2 

+ ( II X || a - l (X,^) 3 - 8} {e^. + ... + e* ± - 1}] 
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+ X k^i " ^ X,e 0^ 6 0i 

k>i 


(X,e K ) e Ni ){t k - (X,e 0 )e 0J , - 


+ 1 ii x r 


- S (x,e ) 3 

j=0 J 


n 

- 6)1 E 
k=0 
kM 


X k e Oi e Ok 


+ 


n 

... + z 

k =0 


X k e Ni e Nlc 


)] > 0 . 


k>i 


We can write this as 
n 

s (v! + B i x i) ^ °» 

where 

A i “ U fc i - (x, e Q )e oi "... - (X>e N )e Ni } 

+ { II X f " " (X,T ) 3 - SHe^ + ... + - 1}] 

j =0 J 

and 

B i = 2[ -o X k Ct i ' <*V*L ' " ( X ’^H )e Ni )lt k ' (X,e O )e Ok 

k>i 

- ... - (X,e N ) eNk } 
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for : 
Case 


+ ( II X || a - s (X,e ) 3 - 6}{ S X fc e e + 

j=0 0 k=0 k 01 0k 

k>i 


+ J 0 ^ e Ni e Hk } ] • 

k>i 

n 

The inequality S (AA® + BA.) > 0 is satisfied if A.X* + BA. > 0, 
i_Q liix ii i i 

i- ® 0, 1, • • » , n « 

1: A i > 0, for some i, 0 < i < n. Assume A n < 0, A h < 0, . . 

A i+1 < 0 and A^ > 0. Then we need to solve 

(1) AjX 3 + B.\ 1 + (A n X 3 + A a _jX 3 _ x + B n _ 1 X n _ 1 + V 2 ^n-2 

+ ®n-2 X n-2 + ••• + A i+i A i+l + B i+l X i+l^ - °‘ 

Thus, we must have 

*1 - U lVn - SWn-l - Vn-l‘,1 - 4 Vh-2 A n-2 


- VnAn-2 - ••• - ViAl - 4A i B i+l A i + l > 

Then choose \ n arbitrarily and X n _ x , X n _ 2 , , X 1+1 as follows: 

sign X = -(sign B.), i+1 < j < n-1. 

In addition, choose X^ to be any solution of (l) and X^ 
to be any solution of A^ 3 + BjX^. = \ fe (A ^ + B k ) > 0, 
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Case 2: If < 0 for all i, choose X^ - 1 and examine 


( 1 ) 


B 3 

n- 


- 4A A n . 
1 n n-1 


If (l) is greater than or equal to zero, we are assured 
of a solution to the equation 

( 2 ) Vl"n-1 + B n-l"n-l + Vn = °* 111611 let "n-1 be either 

B 

solution of (2) and let \ = - — * > J “ 0,1,... ,n-2. 

o A 

If (l) is negative choose = 0, and let = 1 911(1 examine 

(5) ^-2 ' 4A n-l A n-2‘ 

If (3) is positive or zero, then 


V 2 "n-2 + V^n-2 + Vl"n-1 = 0 has a solution - 

B 

Let \ ^ be either of these and let X ^ , 

n-2 J a 

j 

j = 0,l,...,n-3, etc. 


In order for > 0 for some 1 } we must have 

[t. - l (X,e.)e ] 2 + [ || X |P - l (X,e ) 2 - 6] 

1 j=0 J J J=0 

N 

[( Z e *) - 1] >0, or 

j=0 
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„ _ N _ _ N N 

[ II X || 2 - 2 (X,e ) 3 - 6] [1 - 2 e 3 ] < [t . 2 (X,e.)e.J s , or 

J=0 J j=0 Ji 1 .i=o J J 1 


[t i - A 


[ II X II 3 - 2 (x,e ) 3 - 6] < 

j-0 J 


J=0 


But we may write 


1 ' s e ji 3 
0=0 J1 


H 


[t A - 2 (X,e )e ] 3 

N _ 1 j =0 J Ji 

this as 6 > || X (I 3 - £ (X,e.) s + 


6 > E + 


j=0 


[t - 2 (X,e )e ] 3 

j=0 0 


* G ji - 1 
j=0 


N 

E e ,, 3 - 1 
J=0 


Ji 


In the newly computed vector ^ +1 = (X^X^ . . . let \ ± 

be the value of some function at fi ± , i.e cp^C^ \ ± . Then 


E = ifo [X(Pi) " Jo A J^ (Pi) " ^iW p i )P < *• 


Problem ; Determine cp N+1 (3 f ) for some value ^ 0 < i < n, such that 

N+l 

the error obtained by using £ A.cp.(p') to approximate X(p' ) 

0=0 J J 

in the sense of least squares, does not exceed the error obtained 
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N 

by approximating X(p' ) with E A.q>. (£’). 

j=0 J J 

Solution ; Compute A^-^k), k = -1,0, 1, . . . , N, e N+1 and A^ +1 as follows: 


a h+i^ -1 ^ 

w°> 



II %+l ’ ( < %+l ,e j^ e j II 

W -1) V“ 1)( 


a n+ i( n ) “ %+i’V ■ 

N _ 

e N + l = W” 1 * % + i - jS, W J > e j- 

^N+l = ^ ,e N+l^ 

Finally, compute the (N+2)Aj's, j = 0,1, . . . ,N+1, as follows: 

\+l = • a k+i a n+i^ -1 ^ 

= " A H+1 A N+1^ H ^ 

Vi = - W*- 1 * + W-W 1 " 1 * 

+ W n)a n ( n - 1)]} 
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Now let be a such that || p^, - 3' || - min ( || p - p 1 || } . 

° < i < n 

Define the following function 

J 0 Vd (p,) + W(P’)’ 


where 

M(P’) 


h' 


- 2 II Pj. - P' 
L(P it ) 


for 2 || p., - P* || < L^,) 


= 0, otherwise, 

and L^,) = min { || Pj - P ± , || }• 

0 < i < n 
i t i' 

Thus, when p f is chosen, we are able to use the function above to 
approximate X(f3 r ), being assured that the approximation obtained here is 
N 

no worse than the value £ 
j=0 

squares approximating function. 

Writing this multiple of \ as 

§L0i.) - II Pi, - P' II 

■ " “ T r 1 “ > 

\ L(p i( ) 


A^cpj(3* ) obtained by using the initial least 


we see that we have a factor which varies from zero to one as p 1 varies 
from a position on the boundary to a position at the center of the ball 
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{p I ll p 1( - p' ll si } * 


Thus, the factor \ , which was derived in association with the vector & ± , 
is weighted depending on the nearness of p' to * 


Achieser, N. I., Theory of Approximation, Ungar Publishing 
Company, New York, 1956* 

Approximation of Functions, Symposium held at General Motors 
Research Laboratories, Elsevier Publishing Company, 196% 
edited by Henry L. Garabedian. 

On Approximation Theory, Proceedings of the conference held 
in the Mathematical Research Institute at Oberwolfach, Black 
Forest, 1963, Birkhauser Verlag, 1964. 

On Numerical Approximation, Proceedings of a Symposium conducted 

by the Mathematics Research Center, United States Army, at the 
University of Wisconsin, Madison, 1958* edited by Rodolph E. 
Langer. 
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^ General Solution for Minimum Impulse Transfers 
in the Near Vicinity of a Circular Orbit* 

by T. N. Edelbaum 

Analytical Mechanics Associates, Inc. 

Cambridge, Massachusetts 


ABSTRACT 

An analytic solution is obtained for minimum impulse transfer between 
two neighboring low-eccentricity orbits. The orientations of the arguments of 
perigee and the line of nodes are completely arbitrary. The solutions obtained 
are of three different classes. The first class consists of two-impulse nodal 
transfers with both impulses occurring on the line of nodes and having equal 
but opposite radial, circumferential and normal direction cosines for the im- 
pulses. The second class also consists of two-impulse solutions but with equal 
circumferential direction cosines, and equal and opposite radial and normal 
direction cosines. The location of the two impulses for this class is a function 
of the particular transfer. The third class consists of singular solutions with 
a well-defined thrust direction at every point but with an infinite number of 
solutions for the distribution of impulses (or even continuous finite thrusts) all 
having the same fuel consumption. The singular solution can be realized with 
two impulses so that two impulses suffice for all of these optimum transfers. 

^Performed under contract NAS 12-26. 


N67-29373 
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NOMENCLATURE 


a 

a 

c 

C 


e 

e 


x 

e 


y 

H 


x 

i 

y 

K 


p 

R 

s 

T 


u 

X. 

1 

6 

e 

x 

x i 

M 

V 

GO 

O 


Semi-major axis 

Semi-major axis of circular reference orbit 
Cos (6 - e o ) 

Eccentricity 

x component of eccentricity 
y component of eccentricity 
Variational Hamiltonian 
Inclination of orbit planes 
x component of inclination 
y component of inclination 
Gravitational constant 

Direction cosine of the thrust normal to the orbit plane 

Radial direction cosine of the thrust 

Circumferential direction cosine of the thrust 

Magnitude of the primer vector 

Defined by Eq. (50) 

Sin (0 - 6 ) 
v o' 

Defined by Eq. (61) 

Time integral of thrust acceleration 
Defined by Eqs. (1> — (5) 

Defined by Eq. (31) 

Central angle 

Radial component of the primer vector 
Lagrange multipliers 

Circumferential component of the primer vector 
Normal component of the primer vector 
Defined by Eq. (60) 

Defined by Eq. (60) 
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INTRODUCTION 


The theory of optimal space vehicle guidance has much in common with 
classical applied mechanics and celestial mechanics. As might be expected, 
many of the techniques of these classical disciplines have been applied to prob- 
lems of optimal space vehicle guidance. Such techniques include Taylor series 
expansions, linear perturbation theory, higher order perturbation theory, matched 
asymptotic expansions, the method of averaging, and Hamilton -Jacobi theory. 

Many of these techniques have been quite successful for some problems. How- 
ever, optimal guidance problems introduce some special difficulties that are 
not present in typical problems of celestial mechanics. 

One of the usual procedures for transforming an optimal guidance prob- 
lem into a problem in Hamiltonian mechanics is to use the maximum principle, 
or its classical equivalents, to express the control in terms of the adjoint. 

The control is often a highly nonlinear function of the adjoint so that small 
changes in the adjoint may cause large changes in the control. This simple 
fact can cause many of the methods which are successful for problems in celes- 
tial mechanics to break down when applied to problems in optimal guidance. 

This simple fact can even cause trouble with perturbation methods spe- 
cifically developed for optimal guidance, such as the method of "neighboring 
optimal" or "second variation" guidance. This guidance method linearizes the 
state, the adjoint and the control, and can cause difficulty towards the terminal 
point where small changes in the state may require large nonlinear changes in 
the control. The same difficulty can occur with higher order schemes, such 
as n^ 1 variation guidance, because the expansion of the solution in variations 
of various orders may not converge . 
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One method of studying this terminal accuracy problem is to linearize 
the problem about the final state. The resulting trajectory optimization and 
optimal guidance problems may then sometimes be solved analytically. For 
power-limited rockets, the control is a linear function of the adjoint and the 
above difficulties do not arise. Minimum -fuel power- limited solutions have 
been obtained for three-dimensional linearizations about circular orbits (Ref. 
1) and elliptic orbits (Ref. 2). These solutions are for rendezvous in a fixed 
tinle. 


For rockets with constant exhaust velocity, the control is a nonlinear 
function of the adjoint and the problem is far more difficult. Even if no bounds 
are placed on the control magnitude (so that impulses are allowed) and the 
transfer time and terminal positions are left open, only limited (but significant) 
success has been achieved. The coplanar problem for elliptic terminal orbits 
has been partially solved in Refs. 3 and 4. A more complete solution has been 
obtained for three-dimensional transfers in the vicinity of a circular orbit 
(Refs. 4 and 5). The present paper contains additional details on the solution 
of Ref. 4, including the synthesis of the optimal control and the determination 
of the minimum number of impulses required for the singular case. Ref. 5 
represents an independent derivation of the results contained herein and in 
Ref. 4. It includes some consideration of the effects of a slightly eccentric 
reference orbit. 

For a more complete review of the existing literature on optimal orbital 
transfer, the two survey papers prepared under this contract should be con- 
sulted (Refs. 6 and 7). 
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ANALYSIS 


The solution of this problem is conveniently carried out in terms of 
Lagrange's planetary variables . By linearizing the variation of parameter 
equations about a circular reference orbit, the equations of motion (l)-(5) are 
obtained . 


dx 


da/a 


/ a 

_ o / _o 

du du 2 V K l T 


( 1 ) 


dx 


2 

du 


de 

Z 

du 


_o 

K V T 


(2 l T sin 6 - l cos 6) 


( 2 ) 


dx„ 


de 


du 

du 

•*4 

= IV 

du 

du 

dx 

di 

5 

_ X 

du 

du 


— (2 1 cos 0 + sin 0) 

XV I ft 


K S Si “ 9 


a 

K l N COs6 


(3) 


(4) 


(5) 


Following Contensou (Ref. 8) and Breakwell (Ref. 9), the independent variable 
is taken as the time integral of the thrust acceleration. If the thrust is impul- 
sive, this integral is equal to the sum of the absolute magnitudes of the impulses. 
The variable 0 is the angular position in the reference orbit measured from 
the x axis. Both eccentricity and inclination are treated as vectors having 
components along the x and y axes which lie in the plane of the reference 
orbit. The I's are the direction cosines of the thrust in the radial, circum- 
ferential and normal directions . 
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The variational problem is formulated in terms of the Hamiltonian de 
fined by Eq. (6). 

5 dx. 

H = X x i dT 

i=l 


( 6 ) 


Following Lawden (Ref. 10), we shall introduce a primer vector which is the 
adjoint of the velocity vector. The magnitude of this vector will be denoted by 
p and its components in the radial, circumferential and normal directions will 
be denoted by X, fi and V, respectively. The following equations for the op- 
timum thrust direction are derived by means of the maximum principle. 


t 


R 


= A 

P 





v 

P 


(7) 


X 



( - X 2 COS 6 + X 3 Sin 8) 


( 8 ) 





(2X 1 + 2X 2 sin 6 + 2 X 3 cos 0) 


(9) 


V 



( X , sin 0 + X c cos 0) 
4 5 


( 10 ) 


The location of the impulse is given by the value of 8 where p takes 
on its absolute maximum. If there is to be more than one impulse, all of these 
maxima must be equal in magnitude. 


Following Lawden (Ref. 10), Eqs. (8)-(10) will be rewritten in a dif- 


ferent form. 


X 


fir I 2 2 

J i -/ X 2 + X 3 sin < 0 - e o> 


( 11 ) 
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v 


/l 2X 1+ 2 y^r 


Xg cos(e-6 o ) 


( 12 ) 


a 

_o 

K 


XgX^ 




X 2 + X 3 


X 2 X 5 X 2 X 4"*" Xq X- 

2 s “ <®- 0 o ) + 77 2 cos (0- 0 o ) 


X 2 + X 3 


(13) 


tan 0 
o 


(14) 


Equations (11) -(13) are the equations of an ellipse in three -space. This 
ellipse is formed by the intersection of a two-to-one elliptic cylinder parallel 
to the v axis and a plane which passes through the intersection of the cylinder 
axis with the X , p plane. A typical case is illustrated in Fig. 1 which also 
shows the projection of the ellipse on the X , ^ plane. 

There are only three configurations of this elliptical primer locus which 
allow the primer vector to have more than one maximum. The first configura- 
tion is a family of solutions where the center of the ellipse is located at the 
origin (Fig. 2) . In this case the two equal maxima occur on the major axis of 
the ellipse and are separated by 180°. The following equations characterize 
this case, which will be referred to as the nodal case. 


X(e i> = - X(0 2 ) 
)i( e 1 ) = - H(6 2 ) 

v ( 6 ]) = -f(6 2 ) 


(15) 

(16) 


(17) 
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The second configuration, which allows two equal maxima of the primer vector, 
corresponds to cases where the ellipse passes through the p axis and the 
primer vectors again lie in a single plane (Fig. 3) . This case, which will be 
referred to as the nondegenerate case, is characterized by the following equa- 


tions and inequalities 

X 2 X 4 ^3*5 


(18) 

e - e = e - e 1 

2 o o 1 


(19) 

X(6 1 ) = -X(6 2 ) 



= n0 2 ) 

> 

(20) 

I'tfy = -i/(e 2 ) 

( 1 2 2 \ 


4 *lW ^2 + ^3 / 

(21) 

cos (9 X - e o ) - 

2 2 2 2 

CO 

-x 2 x 5 ) - 3(X 2 + X 3 ) 


2 2 2 

3(X 2 +X 3 ) < (X 3 X 4 

2 

- W 

(22) 




2 2 2 2 




(23) 


The third configuration is a combination of the two previous cases where 
the primer locus forms a circle and the primer vector has the same magnitude 
at all points on the orbit (Fig. 4). This singular case is characterized by the 
following equations where the magnitude of p has been taken to be unity. 
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1 

(24) 

*2 X 4 " *3 X 5 

(25) 

"VV ■ <wW- f? 

o 

(26) 

X = ~ sin (6 - 0 ) 
^ o 

(27) 

H = cos (6 - 8 o ) 

(28) 

J^3 

V - Sin (0-6) 

^ 0 

(29) 


The above results for the equal maxima of the primer vector were first 

obtained by purely geometrical reasoning. They have been analytically verified 
by M. Washington in Ref. 4. 

The admissible adjoint solutions have now been determined. The next 
problem is the solution of the two-point boundary value problem to determine 
the optimum transfers corresponding to each adjoint solution. The simplest 
case is that of nodal transfer. The x axis may be aligned with the line of nodes 
between the initial and final orbits. Both impulses must occur on this line of 
nodes because of the 180° central angle separating the impulse locations . 

The total change in the orbital elements is given by Eqs. (30)-(36) using Eq. (17) . 


u 


+ U 2 


(30) 


6 = 


u 


(31) 
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9 

(32) 

(33) 

(34) 

(35) 

(36) 


These equations may be solved simultaneously to determine the total required 
impulse [Eq. (37)] and the magnitude of the first impulse [Eq. (38)] . 



There are two bounds on the region of applicability of this solution. The first 
is given by the requirement that each impulse must point in the positive direc- 
tion of the primer vector . 

(V? * 
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The second inequality is given by the requirement that the primer vector must 
have a maximum at the impulse locations . 

2 2 

Ai * 3Ae y (40) 

The optimal directions of the impulses must still be determined. An 
instructive way to do this is to differentiate the payoff with respect to the state. 

If the maximum magnitude of the primer vector (which is equal to the magnitude 
of the Hamiltonian) is taken as unity, the Lagrange multiplier for each compo- 
nent of the state is the partial derivative of the payoff with respect to that com- 
ponent of the state . 


X 1 “ 0 (41) 
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The direction cosines of the impulses are given by the components of 
the primer vector at the impulses . 



The solution of the two-point boundary value problem for the nondegenerate 
case is more complicated and involves considerable algebra. It is convenient to 
use a different normalization of the Lagrange multipliers. This normalization 
will be retained only to Eq. (67), where the conventional normalization with the 
primer vector equal to unity will be reintroduced. 
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C = cos (8 - 0 q ) 
S = sin (0 - 0 q ) 


( 51 ) 


With these definitions, the direction cosines of the thrust at the impulse loca- 
tions may be written: 


t = A 
R P 


2R 2 S 


/l +R 2 J 4R 2 + ( 1 - 3R 2 ) C 2 


‘TV 


J~L 1 


R C 


74R 2 + (1 - 


2 2 
3R ) C 


V 

" p 


2RS 


yr + r 2 y 4 R 2 Mi-3R 2 )c 2 


(52) 


(53) 


(54) 


The Lagrange multiplier has been eliminated from these equations by use 
of Eq. (21) . The total change in each element of the orbit can now be given by 
Eqs. (55) -(59). 


Ax, 


2 C 


/ 


1 +R 


4R 2 + (1-3R 2 )C 2 


(55) 


Ax rt 


u 


J 1 + R 2 


2(C 2 + R 2 1 

^4R 2 + (1-3 R 2 ) C 2 


(56) 
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u 3 = 2 C S( 1 - 26) 

J 1 +E 2 J 4R 2 + ( 1 -3R 2 ) C 2 

^ x 4 _ 2RS 2 

7l + R 2 ,/ 4R 2 + (X-3R 2 )C 2 

^ x 5 = 2R C S( 1 - 26) 

fl +R 2 74R 2 + <1-3R 2 )C 2 


( 57 ) 


(58) 


(59) 


The x axis has been aligned with 6^ for these equations . This orientation 
will also be abandoned after Eq. (67), when the orientation along the line of 
nodes will be reintroduced. 


The two-point boundary value problem will now be solved in terms of 
the angle between the vectors describing the eccentricity change and the inclina- 
tion change. The following two angles are first introduced. 


tan co 



tan O 


Ax 

Ax, 


(60) 


The variable T defined in Eq. (61) can be determined from Eqs. (56)-(59) . 


(1-26> 2 C 2 S - (C 2 + R 2 )S 

T = tan(io-fi) = i g 

( 1 - 26)( 1 +R )C 


(61) 


Eq. (61) can now be solved for the impulse split in terms of T . 


1-26 = 


o / o 2 2 2 2 

T( 1 +R ) - v T ( 1 +R ) + 4S ( C + R ) 


2SC 


(62) 
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The total change in eccentricity and inclination can now be calculated in terms 
of C, R, S, and T. 



u u 


4(C 2 +R 2 ) +2T 2 (1+R 2 ) -2T7T 2 a+R 2 > 2 +4S 2 (R 2 +C 2 ^ 

2 2 2 
4R (1-3R )C 


(63) 



2 2 2 / o 2 2 2 9 9 

4S + 2T ( 1 +R ) - 2T v T ( 1 +R } + 4S (R + C ) 

2 9 9 < 64 ) 

4R^(1-3R )C 


By multiplying Eq. (63) by R and dividing it by Eq. (55), and dividing Eq. (64) 
by Eq. (55), and subtracting the squares of these two quantities, the following 
result may be obtained. 


C 


2 


2 2 
R (1-R ) 



2 A a 2 2 * 

2R — — - (1+R )(R Ae - Ai ) 

a 

o 


(65) 


This equation allows C and S to be eliminated from the previous equations 
and allows R to be determined in terms of the changes in the orbital elements . 
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The payoff now can be calculated in terms of the orbital elements . 



(67) 


At this point the orientation of the x axis along the line of nodes and the nor 
malization of the Lagrange multipliers with p =1, used in the rest of this re- 
port, are reintroduced. The equation for the total required impulse now be- 
comes Eq. (68). 



( 68 ) 


The Lagrange multipliers are now most easily obtained by differentiating the 
payoff, as was done for the nodal case . 
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The orientation of the primer locus relative to the line of nodes can be found by 
use of Eq. (14). The location of the impulses can be found from Eq. (21) and 
the direction of the impulses from Eqs. (8) -(10). The impulse split between 
the two impulses can be found from Eq. (62) or from one of the boundary con- 
ditions . 
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The payoff for the singular case is independent of the semi-major axis, 
as it is for the nodal case, because the partial derivative of the payoff with re- 
spect to the semi-major axis, \ v is equal to zero. The required value of 
change in the semi-major axis is easily determined by setting X 1 equal to zero 
in the nondegenerate solution for X ^ Eq. (69). 

2 2 2 2 * f a a -2 

Aa = Ae +Ae +~— AiAe -Ai 

x y JJ y 

For changes in the semi-major axis smaller than those given by Eq. (74), the 
solution will be singular (or nodal) and will have the same payoff as the nonde- 
generate solution with the semi-major axis given by Eq. (74) . 


u - 


ft y 7 + *l 


( 74 ) 


(75) 


The orientation of the singular locus may also be found as a special case of the 
nondegenerate case. 


tan 9 


Ae y +7T Ai 

Ae 


(76) 


The optimum thrust directions are given by Eqs. (27)-(29). 


In this singular case, the solution space collapses from five-dimensional 
to three-dimensional and may be visualized in three-space. The simplest way to 
show this is to once again align the x axis with the e o direction of the primer 
locus. The rates of change of the elements are then given by Eqs. (77) -(81) . 


dx 


1 


du 


= 2 C 


(77) 
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<*2 

du 

= f sc 

<*3, 

-2 ^S 2 

du 

" 2 2 S 

<*4 

{K. o 2 
" 2 S 

du 

^5 

JT 

= —sc 

du 


(78) 

(79) 

(80) 
(81) 


It is seen that x„ is linearly dependent on x and that x is linearly 
dependent on Figure 5 is a plot of Eqs. (77), (80) and (81) showing the 
possible changes in the elements as a function of the position of the impulse . 

If only one impulse is allowed, the reachable states lie on the line of intersec- 
tion of a circular cylinder and a parabolic cylinder which form the convex hull 
of the intersection. It is possible to reach points on the convex hull by using 
two impulses (see e.g. Ref. 8). The interesting question is whether points in 
the interior of the volume can also be reached with two impulses . If they can, 
then all minimum impulse transfers between orbits in the near vicinity of a 
circular orbit can be realized with only two impulses . 

A proof that every interior point of the convex hull is reachable with 
two impulses has been suggested by W. D. Hayes and S. H. Lam of Princeton 
University. The geometric interpretation is that a straight line which touches 
the space curve at two points can be passed through every interior point of the 
volume . The proof is constructed by drawing a unit sphere about an arbitrary 
interior point and projecting the space curve onto the sphere. The antipodal 
curve to this curve on the sphere is also drawn. The geometry of the problem 
is such that the curve and its antipodal curve will always intersect. As the 
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intersection is sufficient for a straight line which touches the curve at two points 
to pass through the center of the sphere, two-impulse solutions are always pos- 
sible. Even these two-impulse solutions are not unique, as there are usually at 
least two intersections of the two projected curves. 

The explicit calculation of the two-impulse solutions is difficult and has 

not been accomplished. However, three-impulse solutions are easily calculated 

explicitly. For example, the same value of sin 2 (0 - 0 ) may be used for each 

o 

impulse and the impulse splits adjusted to meet the other boundary conditions. 

Because these solutions are singular, they are also realizable with finite 
thrusts . Any transfer with the optimum thrust directions that meets the boundary 
conditions will realize the minimum fuel consumption. 

The type of solution corresponding to any particular transfer can be found 
from the following diagram. 
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RESULTS 


The total impulse required for any transfer as well as the transfer re- 
gions are illustrated in Figs. 6-11. In each figure, the contours are of constant 
values of Ai/u going from zero at the outer boundary to unity at the origin in 
increments of 0.05. Each figure is drawn for a different value of the angle be- 
tween the eccentricity change vector and the inclination change vector. Figure 

6 corresponds to co-axial transfers where the solution is known to consist of 
inclined Hohmann transfers. There is no singular solution in this case. One 
set of Hohmann transfers is a special case of the nodal transfers while the other 
is a special case of the nondegenerate transfers. This becomes clearer in Fig. 

7 where the singular solution appears near the outer boundary of the figure. The 
nodal case is the triangular region adjoining the singular region. The nonde- 
generate case occupies the rest of the figure. 

As the angle between the eccentricity change and the inclination change 
increases, the nodal region rapidly decreases in size while the singular region 
grows. Finally, in Fig. 11 the nodal region has completely disappeared. 
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SUGGESTIONS FOR FURTHER RESEARCH 

The solutions obtained herein are quite general and contain a number of 
new results. They naturally suggest a number of possibilities for further re- 
search, particularly for nonlinear transfers . The following paragraphs will 
briefly describe some of these research areas. 

1 . The existence of more general nodal transfers than the Hohmann 
transfer in this problem suggests a general investigation of nodal 
transfers . 

2. The singular solution obtained herein may be the limit point of a 
nonlinear singular solution like the Lawden spiral. Unlike the 
Lawden spiral, half of this curve satisfies the necessary condi- 
tion derived by Kelley, Kopp and Moyer (Ref. 11) and by Rob- 
bins (Ref. 12) . 

3. Either or both of the infinitesimal impulses of the present solu- 
tions may be allowed to become finite. The resulting nonlinear 
transfers may be investigated by the methods of Moyer (Ref. 13) 
or Winn (Ref. 14) . 

4. The consideration of higher order terms will remove the de- 
generacy of the singular solution and will introduce three -impulse 
solutions (if not even more complex ones) . The next step is the 
consideration of quadratic terms. 

5. The terminal guidance problem can be considered since the sen- 
sitivity of the transfer to various errors can be analytically de- 
termined. Finite thrust guidance might be approached as Rob- 
bins does in Ref . 15. 
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Rejection to Infinity in the Problem of Three 
Bodies When the Total Energy Is Negative* 


by D. C. Lewis 

Control Research Associates 
Baltimore, Maryland 


1, Introduction 

In the problem of three bodies, we take the kinetic energy 
relative to the center of mass of the three bodies and we take the 
zero level of potential energy as occurring when the three bodies are 
infinitely far from each other. This makes the potential energy always 
negative for any actual configuration, and the total energy must then 
be negative also provided that the velocities of the three bodies relative 
to the mass-center are sufficiently small for the given configuration. 

In this context we consider motions for which the total energy has a 
preassigned negative value - K(K > 0) and for which the length of the 
angular momentum vector (taken relative to the mass center) also has 
a preassigned value f > 0. 

^Performed under contract NAS 12-93. 
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A useful quantity for measuring the dispersion of the three 
bodies is the so-called Lagrangian inertial radius R relative to the 
mass center. Its definition will be deferred to the next section. We 
merely note here that R becomes infinitely large if and only if the 
perimeter of the triangle formed by the three masses becomes infinitely 
large, and R is small if and only if the perimeter is small. 

G. D. Birkhoff has shown that if R is at any time t Q less 
than a sufficiently small positive quantity Ry, then one of the three 
bodies will recede infinitely far from the other two as (t - t Q )+ <®, 
while the distance between the latter two will remain bounded. 

It is clear from his work that Uy depends on the three masses 
niy, m^, m 2 as well as on the assigned values of K and f, so that 

R 0 = V K ' f ' V V V * 

But he did not give an explicit estimate of Ry. The main purpose of this 
paper is to calculate such an estimate. The final result is formulated 
with some precision in Theorem 10. 

We also wish to clear up some obscurities, if not actual errors, 
in Birkhoff's treatment. For instance, Birkhoff never proved the statement 
about upward concavity of the part of the curve R = R(t) for which 
R < f/(2 1 ^ 2 K 1 ^ 2 ), cf. Reference [1] p. 278, l. 1. He only proved that 
the second derivative was positive at points where the curve was horizontal. 
For such reasons we have discussed this theorem in much greater detail 
than was done by Birkhoff. Our version of this part of his work appears 
as Theorem 3. 
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We have ignored the possibility of collisions in this paper. 
Since f > 0, the case of a simultaneous collision of all three bodies 
*is ruled out by the conjecture of Weierstrass, which was rigorously 
proved by Sundman, The remaining difficulties afforded by the possible 
occurrences of binary collisions are readily eliminated with the help 
of a regularizing parametrization due also to Sundman. 


2. Notational Prelude. 


The nine second order differential equations for the three 


body problem may be written in the form of three vector equations 
d 2 

(2.1) n ' ".y 1 * [ “ » i = 0, 1, 2 f 


1 dt 2 


dU 

3q . * 


where q^ is the vector with components x^ # y^, z^. The mass nu 
is thus regarded as having rectangular coordinates x^, y^ t i ^ in an 
inertial frame of reference. The (scalar) force function is 


m l m 2 


m 0 m l 


where r^ is the distance between the masses nu and Here (i, j , k) ^ 

is a permutation of (0, 1, 2), and ^ • |q^ - qj * [(Xj-x k ) 2 +(y..-y k ) 2 +(i..-z k ) 2 ]^ 

The equations (2.1) are well known to admit a total of ten 
elementary first integrals, the first six of which express the conservation 
of linear momentum and permit us to use a frame of reference whose origin 
is at the center of gravity of the three bodies. We may thus assume that 


( 2 . 2 ) 


m Q q 0 * m l q l + m 2 q 2 2 ° * 
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The energy integral is written in the form 

2 7 \ IqJ 2 - U - k , 

1 

where the dot is used to denote differentiation with respect to t and 
where the integration constant K is the negative of the total energy. 

The three integrals of angular momentum appear as a single 
vector equation 

2 m i (q i X qp = c . 
i 

The length of the vector c of total angular momentum will be denoted 
in the sequel by f. 

The Lagrangian inertial radius R may be defined by the 

equation 

R 2 = Z m. |q.| 2 , 

i 

but, if (2.2) is assumed, it is not hard to see that we also have 
R 2 = M~ 1 (m^m^r 2 + m^r 2 + m^r 2 ) 9 

where M, as in (2,12), is the sum of the masses. 

Lagrange’s well known identity is to the effect that 
2 2 

(2.3) ^4- = 2(U - 2K) 

dt 

We refer to the following as Sundman’s inequality 

(2.4) R 2 ♦ 2R R* ♦ 2K > f 2 R -2 
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which is well known. The auxiliary function of Sundman is defined by 
the equation 

(2.5) II = R R 2 + 2 KR + f 2 R _1 . 

Its derivative with respect to t is seen to be of the form F R, where 
F = + 2R R + 2K - f** R ^ ^ 0 because of (2,4), We thus note the 

important fact that, on any interval for t on which R is monotonic, 

II is also monotonic in the same sense. 

The formulas thus far used have the advantage of being 
symmetric in m^, m^ ( and q^, q^, In other words they are 

invariant under any permutation of the subscripts (0, 1, 2), On the 
other hand they have the disadvantage of not reflecting the full potentiality 
of (2,2) for reducing the number of unknowns. One way of accomplishing 
this is due to Lagrange, but a single reduction of this sort necessarily 
sacrifices the desirable symmetry noted above. Consequently we contemplate 
the whole class of such reductions in the following way. 

Let (i, j, k) be any permutation of (0, 1, 2), Let the 

vector q, with components (x, y, 2 ), determine the position of the 

mass m. relative to m. . so that 
J 

(2.6) q “ q j ’ q i * 

The center of gravity of nu and iik , relative to the origin of the 

original frame of reference, is evidently at the point corresponding to 

the vector a.q. + a.q., where a. *= m. (m. ♦ in.)”* and a. = m.(m. ♦ m.)“* . 

IX * 1 i v 1 y j J 1 J 
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We now use a vector s with components (£, n, O to determine the position 
of the third body m k relative to the Center of gravity of and m y 

so that 


(2.7) s * V*i * a j q j 

Assuming that the origin of the original frame is at the center of 
gravity of all three bodies, 50 that m.q^ ♦ 4 m k q k * 0* which is 

simply another way of writing (2.2), we may easily find q^, q^, and 
q k in terms of q and s. In fact we have 

q i “ * a j q * \ M ~ 1 s 

(2*8) fjj * * “ \ M 1 5 

q k * (m i + mp M” 1 s 

where M, as previously, denotes the sum Of the three masses. 


If we set 


( 2 . 0 ) 


m.m . 

*JL 


m. + m. 
1 J 


and y 


(m i ♦ 11^) 

M “ 


it may be shown that the equations (2.1) are equivalent, with due regard 
to (2,2), to the equations 

.. 9U «• 3U 

m * * W S “ as ’ 

The energy integral may be written in the form 
(2.10) j m|q| 2 ♦ ip| s| 2 » U - K 
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and the Ugrangian inertial radius R, previously introduced, satisfies 
the equation 

(2.11) R 2 • m x 2 * u p 2 , 

where r ■ |q| » / x 2 * y 2 + z 2 and p » |s| » / £ 2 ♦ n 2 + ? 2 . 

Most of the above formulas, beginning with (2,6), depend 
heavily upon the particular choice of the permutation (i, j, k) of 
(0, 1, 2). It will be essential to have upper and lower bounds for both 
m and u which are independent of this permutation. We will, in fact, 
calculate such bounds in terms of the following four symmetric functions 
pf m^, i and ^ 2 * 

M * m^ + 

P » m 0 m 1 m 2 

( 2 . 12 ) 

* 

m » minimum of the three masses . 

m m greater of the two smallest masses , 

2 .2 

From (2.9) we find that dm/ dm. * m.(m. ♦ m.) > 0. Hence m is an 

* 3 • J 

increasing function of when nu is held fixed. On the interval 
"ij < ^ < * ») it takes on its minimum at nu » , at which point 

* ■ i ■ j i i m ; and as + *> , m tends to its least upper bound 
m j , and, since under present conditions £ kk, ru cannot exceed 
the greater of the two smallest masses, which we denote by m. We 
conclude therefore that 
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(2.13) 


■j m ^ m < m 


if m. > in.. But. since m is symmetric in m. and m. (and so 

i - j * i j 

★ 

are m and m) , we see that (2.13) holds also even if nu < m_. . It 
is seen from (2.9) and (2.12) that y = P m * ff*. Hence it follows 
at once from (2,13), that 

(2.14) < y <, y , 

where 

- 2P 


(2.15) 


y = — ■ and y = 

m M m M 


3, Fundamental Results . 

Theorem 1 . In case K > 0, the least of the three mutual distances, 

2 

Tg, r^ , r^ can not exceed M /(3k). 

Proof. From the energy integral it is clear that 


K < U 


2 2 

where r = min (r Q , r^, r^) . But M = (m^ ♦ ♦ m^) 

12 2 12 2 12 2 
3 yO 0 + n> x ) ♦ j(m Q + m 2 ) + j(m 1 + m 2 ) + 2(m Q m 1 + m Q m 2 + m 1 m 2 ) 

> 3(1^1^ + m^m 2 + m m 2 ) . Thus we have 


K < 


m 0 m i + m 0 n 2 + m l m 2 ^ it 
— Y =» 3r 


so that 


as desired. 


3K 
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We continue throughout the rest of the paper to assume that 
K > 0. Moreover, as in the proof of the next theorem, we shall always 
use r to denote the least of the three mutual distances and we shall 
use p to denote the distance from the center of gravity of the two 
mutually closest bodies to the third body. This notation is consistent 
with the notation of (2.9) and (2.11), since we merely have to take rru 
and jik as the two closest bodies. 


Theorem 2. If 


:>r /= 

R > ^ / m + 4 y 


1 2 - 1 

then p > ■=* M K . 


Proof . Evidently, by (2.13) and (2.14), 


M 

R ^ / m + 4 y , 


Since R =mr + y p , we have, from Theorem 1, 


2 _ 1 „2 in _2 M 4 . , m M 4 4 M 4 


p - - R - - r > ~ 

Vi U — n ,-2 


(n + 4 w ) - - 7 =■ 

9 K‘y * 9 K 9 IT 


so that 


P > 


2 M 

3 K 


as stated. 


In the sequel it is convenient to refer several times to the 
following elementary lemmas. Although they are essentially well known, 
it is more convenient to give the proofs here than to refer to the 
pertinent literature. 
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Lemma 1 , Suppose that F(t) is differentiable for t > a and that 

lim F (t) « b, where a and b are constants* Then there exists a 
t**» 

sequence t. , t n , ... such that lim t » • and lim F'(t ) « 0. 
n 1 2 n+<*> 11 n-*® n 

Proof . Assuming n > a, we may define t^ # by the mean value theorem, 
to be such that n < t^ < n ♦ 1 and F*(t n ) 3 F(n ♦ 1) - F(n), and 
the lemma follows at once since lim [F(n + 1) - F(n)] * b - b * 0. 

nr*® 

Lerama 2 . Let F(t # P) be a continuous real valued function of its two 
real arguments t, P; and let it satisfy a local Lipschitz condition in 
P for t Q <, t and P > 0, Let R(t) be a positive differentiable 
function of t, defined for t £ t^, such that 

(I) R'(t) < F[t f R(t)] 

Let Q(t) be positive and satisfy the differential equation 
Q'(t) * F[t , Q(t) ] 

on the interval < t < t^ and suppose also that Q(t Q ) * P(tp) . Then 
q(t) >_ R(t) for t Q < t < t x . 

Proof. Let t^ be aj wy real number between t Q and . Then by the 
existence and continuity theorems for differential equations, we may define 
a function P(t) * P(t, e) which satisfies the differential equation 

(ID P'(t) * Fjt, P(t)] ♦ c , 

and the initial condition 

PCV * R(t 0 ) . 
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this definition being valid for t Q < t < t 2 , at least, if | e | is 
sufficiently small, depending perhaps on the choice of t 2 * 

Moreover 

(HI) lim P(t, e ) * Q(t) 

e-*0 

on the interval t Q <* t £ t 2 . 

We first prove that, if e is positive, then 

CIV) P(t, e) > R(t) for t Q < t < t 2 , 

Suppose that this were not so. Then there would exist a point t*, such 
that t Q < t £ t 2 , where the function w(t) * P(t, e ) - R(t) is not 
positive. Also, since w(t 0 ) = P(t Q , e) - R(t Q ) = 0 and 

CV) W(t 0 ) = P*(t 0 , c) - R'(t 0 ) > [F(t 0 , P(t 0 )) ♦ E ] - F(t 0 , R(t 0 )) * c > 0 

** ★ 

there is a point t > t Q , but not as great as t , such that w(t) > 0 
on the interval tQ < t ^ t • For otherwise we would have w(t) non- 
positive at an infinite number of points in any neighborhood of t Q , 

necessitating w'(t 0 ) £ 0, contrary to (V). Thus, since w(t**) > 0 

* 

and w(t ) < 0, there is a non-empty set S of points on the interval 

** 

t < t < t 2 where w vanishes. Since w is continuous, S is closed. 

Let t ■ the greatest lower bound of S. Then t Q < t**<7_< t 2 and 
w(t) « 0. That is, 

(VI) P(7, e) * R(t) , 

But w(t) > 0 for t^ < t < t , Since w is differentiable at 7, it 
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follows that w'(t) 4 0, which contracts the fact, following from 
(VI), that 

W(t) = P'(t,€) - R f (t) > [F(t, P(t)) ♦ e ] - F(t, R(t)) = e > 0 . 
This finishes the proof of (IV) , 

By letting e -► 0 in (IV) , we find from (III), that 
Q(t) > R(t) for < t < t^. But, since t^ is arbitrar>’ on the 

open interval from t^ to t^ and since both Q(t) and R(t) are 
continuous, the above inequality remains valid on the whole closed 
interval from t^ to t^, as we wished to prove. 


Theorem 3 . In case f > 0, K > 0, any connected part C of the curve 

- 1/2 

R = R(t), for which R < f(2K) , consists of an arc with a single 

minimum. If R = R at this minimum and if 0 is any number between 

2 

0 and 1, the curve rises on either side until R > 0 f /(2K R Q ) with 
slope R = (d R/dt) at least as great in absolute value as demanded by 
the inequality 


R “ ^0 

V 2K) 


R 


at every intermediate stage. It is hereby implied that R^ > 0 . 


Proof , If *R*(t) ^ 0 and R(t) = 0, the Sundman inequality, (2.4) , 

R 2 + 2R k + 2K £ f 2 R“ 2 , shows that R(t) > i f(2K) *^ 2 . We thus see 

1/2 

that, as long as we confine attention to C (for which R(t) < f(21C)“ ) , 

we must have R > 0 at all points where R = 0, This means in particular 
that R cannot be a constant on C or any sub arc of C. It also means 
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that R can vanish at not more than one point of C. We proceed to 
show that such a point actually exists. 

Otherwise, since ft(t) is continuous and never zero, R(t) 
would be strictly monotonic (in one sense or the other) on C, and 
hence the Sundman function 

H = R R 2 + 2K R + f 2 R -1 


must be monotonic in the same sense. This shows that R cannot decrease 
monotonically to 0 for either increasing or decreasing t, tending to 
either a finite or an infinite limit, for this would make II tend to 
♦ oo, so that H could not possibly be monotonic in the same sense as R. 
Hence let R Q > 0 be the greatest lower bound of R(t) on C. Then, 
since R(t) is presently assumed to be strictly monotonic on C, we 
must have either 


(3.1) 
or else 


lin R(t) * R q > 0 
t-*» 


(3.1*) liro R(t) = R > 0 , 

t+-* U 

and since the differential equations are invariant under change of sign 
of t, it will be sufficient to confine attention to (3.1). We wish to 
show that (3.1) leads to a contradiction. By Lemma 1 there is an infinite 
sequence t^, t 2 , ... tending to «, such that lim R( = 0. Hence, 
from (2.5), we have ^im II (t n ) = 2K ♦ f 2 R ^ . But, since II is monotonic 

in the same sense as R and is never negative, we know that lim H(t) exists. 


85 



TRAJECTORY ANALYSIS AND GUIDANCE THEORY 


It follows therefore that lim Ii(t) * 2K R^ + f 2 R^*. Moreover 
2 - 1 

|l(t) > 2K R 0 + f R^ for every finite t, since this corresponds to 
the monotonic sense in which H tends to its limit. From the definition 
(2,5) of H(t), we find therefore that 

R R 2 * 2K R + f 2 R' 1 > 2K R Q ♦ f 2 R* 1 
which is equivalent to 

•2 R * ,( 0 f 2 

(3.2) K V 2K) • 

And, since R(t) is presently assumed to decrease monotonically to R Q 
as t « # we also have 

(3.3) R(t) < 0 
on C, 


We wish to show that (3,1), (3.2), and (3.3) are mutually 

* 

inconsist ant. To do this we take any fixed t corresponding to a point 

* * 

(t , R(t )) on C and define a function Q(t) by means of the differential 


equation, 


(3.4) 


Q(t) * - [• 


- V? f f 2 


Q(t) ] l R 0 Q(t) 

* * 

and the initial condition Q(t ) ■ R(t ). Such a definition can be effected 
by the usual existence theorems for differential equations at least for 


- 2K] 


some sufficiently short interval t £ t < t « It follows from (3.2), 

(3,3), and (3.4) that 0 > Q(t) £ R(t) whenever Q(t) * R(t) > R Q , in 


particular when t *« t , Hence by Lemma 2 the two curves can never cross 
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for t > t and moreover we must have Q(t) £ k(t) > R Q# Since R(t) 

is presently assumed to be defined for t* < t < «, and since, from (3,4), 

$(*) 4 it is clear that lim Q(t) • Q**, say, must 

t+t** 

exist ahd must exceed R Q . Furthermore Q* # <, Q(t*) « R(t # ) < f(2K)“*^ , 

since (t\ R(t*)) is on C. Thus R Q < Q** < f(2K)' 1/2 , so that both 

Q - and f 2 /R q Q - 2K are positive when Q « Q**. it follows by 

repeated application of the existence theorems for differential equations 

that the definition of Q(t) can be extended over the whole infinite 
* 

interval t <t < ■, and that over this whole interval it is monotonic 
non-increasing and bounded from below by R . Hence 


(3.5) q o . lim Q(t) 

exists. We also see from Lemma 1 that there exists an infinite sequence 

of points t., t ? , ... tending to ♦ such that lim 6ft ) » 0 . 

* * n->** v n' 

But from (3.4) and (3.S) we know that lim Q(t) exists (even when we do 

t-K® 

not restrict ourselves to such a sequence). We are thus justified in 
writing 

(3.6) lim Q(t) » - 2.]T . 2K] T - 0 . 

^0 0^0 

Since the point (t , Rft )) is on C and since both Q q and R Q do 

* 

not exceed R(t ), as follows from the monotonicity of Q(t) and R(t) ? 
we have 



f 2 _ 

R(tV 


> 2K . 
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From this result and (5.6), it follows that Q 0 = R Q . We have thus 
established for the function Q(t) that 


(3.7) 

and that 

(3.8) 


lim Q(t) = R q 


lim Q(t) = 0 . 

t “H*> 


that 


(3.9) 


An elementary calculation based on (3.4), (3.7), and (3.8) shows 


• * 1 f ^ if 2 

lim Q (t) - 2T (*2 - 2K) > ( - 2K ) > 0 

2R 0 R 2 0 R(t ) 


But from (3.8) and Lemma 1 there exists a sequence of points tending to 
60 such that Q(t) tends to 0. This contradicts (3.9), and hence we 
have finally shown the inconsistency of (3.1), (3.2), and (3.3), 

This finishes the proof of the fact that C contains just one 
point (t , R(t (} )) where R takes on a minimum value R Q = R(t^) . It 
has also been established that R^ > 0. 

Evidently the curve R = R(t) must rise on either side of the 
minimum point (t^, Rq) on ^ e ^ ier indefinitely or until it reaches 
a relative maximum point (t^, R^) , where, of course, R(t^) * 0, 

R(t x ) 4 0, so that this point (t p R^) is beyond the upper bound for 
points on C. Thus R > f(2K)" 1//2 > R Q . From the Sundman inequality, 
H(R(t Q ), R(t Q )) ^ H(R(t 1 ), R(t x )) together with the fact that R(t Q ) 
and R(tp are both zero, we find that 
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(3.10) 2K R q ♦ |-< 2K R x ♦ (Rj - Rttj)) . 

This can be written in the form 2K R 2 R y + f 2 Ry - 2K R 2 R^ - f 2 R^ ^ 0 
or 2K R x R 0 (Rj - R 0 ) - f 2 (Rj - R Q ) >, 0. Since (R - R q ) > 0, we 
find that 2K R^ R y > r . In other words Rj >_ f 2 /2K R Q . Thus, R 
must rise ter at least this value, if there is a relative maximum. 

If there is no relative maximum and if R still does not 
2 -1 

attain the value f (2K R Q ) , it is clear that R(t) is monotonic and 
2 - 1 

bounded by f (2K R ) , so that it tends to a limit as t -► ± «. This 

2 - 1 

limit is greater than R Q but not greater than f (2K R Q ) , So we 
may write 

C3.ll) R q < lim R(t) = R ^ f 2 (2K R^* 1 . 

t-*» 

We hereby restrict attention to the case t -* «►, as we are entitled to 
do because of the invariance of (2,1) under change of sign of t. Let 
R^ be any number such that R^ < < R, and let 6 be any number between 

0 and 1; then set 

I I - I 

(3.12) e ■ f (R 2 - R 0 ) 2 (1 - 6) 2 R q 2 R _1 > 0 . 

From (3,11) and Lemma 1 we can find t^ > t^ such that 0 <* R(t^) < e 
and such that <* R, where R^ = R(t^). Using the Sundman 

inequality as before we now get a slightly modified form of (3.10), namely 
R^ e 2 ♦ 2K + f 2 >, 2K R Q + f 2 R Q This leads to the inequality 
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„ r K i E > r* C 

R 1 - 2K K 0 ' 2K(U 1 - 2KR Q 2K(P > 2 - R^) 

Using the expression for e given in (3.12), we find that 
R ! L f 2 6/2K R q , as we desired to prove. 

For any t between t^ and t^, either in the case just 
treated or in the previous case where R(tp .> f^( 2 K U Q ) we of course 
have R(t) > U Q ; and hence, the Sundman inequality, with the monotonicity 
of R(t) , leads to the result that R(t) R(t ) 2 + 2KR(t) * f^R(t) 1 ^ 2KR Q ♦ f R 
From this, we get after an elementary calculation the inequality given in 
the statement of the theorem, the proof of which is now complete. 


Theorem 4 . If Ry < nin[f( 2 K)' 1/2 , 3f 2 e 2 _1 tf 2 (m + 4T)‘ 1/2 ], then p 

2 

attains a value at least as great as 2M / (3K) . 

Proof . By Theorem 3, R attains a value > ef 2 /(2K R Q ) which cannot be 
less than 


9 f* 




2 k [. 


\jf _ 

(« + 4u) 


1 _ 

2 


Hence, by Theorem 2, the corresponding value of 


p _> 2M 2 /(3K). 


Theorem 5 . Let X = 2M 2 / [ 3 (K m*) 1/2 ], then |r r| <, X . 

Proof . From r 2 » x 2 ♦ y 2 * z 2 and r * x(£) + y(£) + z(^) we find 
from the Cauchy-Schwarz inequality that r 2 4 x 2 ♦ y 2 ♦ z 2 , which by the 
energy integral (2.10) is less than 2m ^ U, Since r is the least of 
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. the three mutual distances, we have 
« 2 ^ (m^m^ + m Q m 2 ^ m ^ in 2^ 


As shown in the proof of Theorem 1 , + m 0 m 2 + < in 2 . Hence 


r 2 j.2 < 2 M 2 r _ 4 M 2 r 


3m 


3 m 


by (2.13). Now, by Theorem 1, r _< M 2 (3K)“ X . Hence 
,,4 


r 2 r 2 < 4JQ 


9 m K 


so that 


2 M 


3(m K) 


1/2 


X, 


as desired. 


Theorem 6 . If R R > B(R) 1/2 , where B(R) = (4 M 1/2 m 3/4 R 1/2 + X m) 2 , 
then p p > 4 M 1 ^ 2 p 1 ^ 2 . 

Proof . By (2.11) p p 2 < R 2 , so that 

i i I 

(3.13) R 2 > p 4 p 2 . 

Also by differentiation of (2.11), we have ppp - R R - m r r. Applying 
Theorem 5, we see that - m r r ^ - xm, while 

R ft > 4 M 1/2 y 3/4 R 1 / 2 ♦ xm > 4 M 1/2 y 3/4 R 1/2 ♦ A m , in accordance with 
(2.13) and (2.14). Hence ypp > 4 M^ 2 y 3 ^ 4 R 1 ^ 2 . Referring back to (3.13), 
we see that ypp > 4 M 1 ^ 2 y p 1 ^ 2 , whence the stated result follows on 
division by y. 
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Theorem 7 . Let 


and let 


If 


e - min I fC^) » 


ef 


3 0 f 

“ .,2 , — 1/2 

2 M (m + 4 ii) 


2 K e 


Rg min 


o_ 

c » o * 


f 2 o 


2 * 2[3(o) ♦ K a 2 ] 


then R takes on the value a at some later time t = and the 

* . , 2 *2 

corresponding value of R is such that o R > B(o) « 


Proof. If 


R o <*<- — < f 


“ m /2K 


we know from Theorem 3 that R will attain a value 

f 2 / 

1/2 


£ 0 f 2 / (2K R ) >, 0 f 2 /(2Ke) = a . R also takes on the value 


R < e < 


e 


SIX 


say. But then 


e r 

2K 


so that 


. 6 f 2 Of 2 
n ' WJ- 2KT ‘ ° * 

Hence R takes on values both less than a and greater than (or equal 
to) a . Hence, since R(t) is continuous, it takes on the value a at 
some t * t D . 
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By the second part of Theorem 3, we see that at t = t 

. V 

where R * a, we must have 


a 2 R 2 > Co - R 0 ) (f- - 2Ko) > | 


as we wished to prove. 


Y- 


P[B(ir) 


” K^]) 


- 2 K 


B(<y) 


Theorem 8 . Under the hypotheses of Theorem 7, the value of p at 
is > | |r- , while that of R is > M^SK)" 1 (m + 4 y) 1/2 . 


t « t 


«0 


Proof . 2 

0 = frr and e < — 2 - ! ef _ m • Hence 

E 2M 2 (m + 4 7) 1/2 


i 



ef 2 

_ 'A 

P(m 4 y)^ 

2K 

3 0 f 2 '1 


3K 

2M 2 (in + 4 7) 1/2 J 




Moreover R = a when t = t Rfl . Hence at t = we have 

r > m 2 (" 1 * y) 1/2 
3K 

2 

Hence, by Theorem 2, we have p > i when t » to 

— 3 K Kg • 

2 

Theorem 9. If 0 > 1 1_ f then p’ > - 8 M p" 2 . if for any such value 

of p, we have ^ > 4 H 1/2 p' 1/2 , p win constantly increase without 
bound. 

IVe omit the proof of Theorem 9 because the proof in [1] is 
entirely satisfactory. 
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Theorem 10 , In the three body problem with masses m^, m^, suppose 

the total energy - K to be negative; so that K > 0. Suppose also that 
the angular momentum c of the system about its center of mass is not 
zero; so that f = |c| > 0. Let R be defined by the formula 



where M * m 0 + m i + m 2 ^ where r i is the distance between m j 
and n^, j t i, k i i, j i k. 

I, Then the minimum of these three mutual distances is under 

all circumstances not greater than M (3K) and there exists a positive 

number R Q such that, if R <, at any time t - t Q , there exists 

— 1/2 

a later time t^ such that U(t^) > ♦ * u) where m is the 

— _i * * 

greater of the two smaller masses and y * 2 m^m^ M /m , m being 
the least of the three masses. 


II. Por t >, t^, the pair of masses closest together retain 

their identities and Ji® R(t) = so that two of the r. become 

2 - 1 

infinite while the third one remains bounded (by M (3K) ) as t -► «. 


C 4 M 1/2 " 3/4 R 1/2 ♦ Xm) 2 


III. A number R 0 for which the above is true may be calculated 
2 * 1/2 

as follows: Let X = 2 M /[3(K m ) ] and 

Choose any positive number 6 < 1. Let 

r2 


and let 


= min 
0 f 2 

TXT 


l 

0 J 




3 e r 


2 M 2 (m ♦ 4 u)"^ 


Then finally take 
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R q 3 rain e 


o_ f 2 a 

2 ’ 2[B(o) + K o 2 ] 


Proof. The assertion under I to the effect that the least of the 

three mutual distances never exceeds M (3K)’ A follows from Theorem 1. 

The other result under I follows from Theorem 8, assuming that e and 

R y are determined as described under III, so that the hypotheses of 

Theorems 7 and 8 are both fulfilled. The t^ of the present theorem 

can, of course, be taken as the t D of Theorems 7 and 8. 

R 0 

From Theorem 7 we also have 


(3.14) R 2 R 2 £ B(R) 
when t = tj, since R(tp = o . 

Let nu and nu be the two bodies closest together at t = t^ 
and let their distance apart at any time t be r = r(t) , Let the third 
body be distant p (t) from the center of gravity of nu and nu . 

Then from Theorem 8 wc know that at t = t^ we must have 

(3.15) p > 2 M 2 /(3K) ^ 2r . 

Now if, for t > t x , p never decreases, the inequalities (3,15) persist, 

since, as long as r(t) remains the least mutual distance, 2r can never 
2 

exceed 2M /(3K) whereas both r^ and r^ exceed p-r>2r-r=r, 
and for r^ (for instance) ever to become less than r, it would be 
necessary (by continuity) for and r to be equal at some t* > t 9 

But this is impossible, since the posibility of changing the first inequality 
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sign in (3.15) to an equality sign is ruled out. This means that the 

pair of masses closest together retain their identities for t > t^ , ’• 

at least if p never decreases. 

That p does indeed never decrease for t > t^ but actually 
increases without limit may be seen as follows: From (3,14) and 

Theorem 6 we see that at time t = we must have 

I - I 

(3.16) p > 4 V? p ^ 

We now combine (3.15), (3.16) and Theorem 9 to show that 

(3.17) lim p(t) = “ . 

t-H*> 

Finally, the proof of Theorem 10 is completed by the remark that the 
result 


lim R(t) * » 

t-Ko 

follows from (2,12), (3.17), and Theorem 1. 


Reference 


[1] George D. Birkhoff. Dynamical Systems, American Mathematical 
Society Colloquium Publications, vol, IX, Chapter IX. 
Especially pp. 275-282. 
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An Algorithm To Obtain Series Expansions 
in the Three-Body Problem 


by P. Sconzo 

/BM, Cambridge Advanced Space Systems 
Cambridge , Mossac/tvse/ts 


SUMMARY 


N67-29375 


This report on Contract NAS 12-98, "Modern Celestial Mechanics, " 
awarded by NASA-ERC to IBM- Cambridge Advanced Space Systems, deals 
with the analytical solution of the three -body problem by means of convergent 
power series expansions. The method followed uses only simple algebraic 
tools and fully exploits the properties of Levi-Civita's regularizing transfor- 
mation. Stimulated by the suggestion made some years ago by Vernic (1955) 
to take advantage of this transformation for purposes of practical computation, 
we present here, for the first time in the long history of the three -body 
problem, a workable algorithm for constructing recursively its power series 
solution in terms of Levi -Ci vita's regularizing variable. This method solves 
the three -body problem formulated in its utmost generality, since no restric- 
tions at all are made on the order of magnitude of masses and distances and 
none of the three bodies is restricted to move along a conic section orbit. 
Besides, the reference system used is a tridimensional inertial one. 


After an introductory section dedicated to the historical background of 
the problem, the second section deals with the discussion of the algebra in- 
volved in the derivation of the time series solution. In a third section, the 
radius of convergence of this solution in the neighborhood of a non-collision 
point is determined according to Sundman. Section number four describes 
the property of regularizing variables with particular emphasis on Levi- 
Civita's variable u. Finally, in the fifth and last section the power series 
solution in the new variable u is constructed by a procedure of successive 
approximations in which only elementary algebraic operations are to be 
performed on polynomials. 


In the conclusion, a comment is made about the merits of the algorithm 
described in Section V, and the items now being developed to implement this 
investigation are listed. 
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I. INTRODUCTORY BACKGROUND 

The first remarkable contribution to the analytical solution of the 
general problem of the motion of three bodies subject to mutual Newtonian 
attractions was made by Painleve (1896). Here, the designation ''general 
problem" is used in contrast to "restricted problem, " which characterizes 
a special case of the three -body problem, namely, that when one of the three 
masses is negligible and the two bodies of finite mass revolve around one 
another in circular orbits. 

In various papers and in his masterly lectures at Stockholm ( 1 897) , 
Painleve 1 demonstrated that: a) the motion is regular and the coordinates 
may be represented by series of polynomials in t (time) in any interval in 
which no collisions occur; and b) the motion ceases to be regular in either 
of the following cases, when the three bodies collide at the same point or 
when only two of them collide, their distance from the third body remaining 
a finite quantity. The analytical characterization of the occurrence of col- 
lisions was later examined by Levi-Civita^ (1903.-04) and his disciple 
Bisconcini^( 1 904K They also investigated qualitatively the holomorphic re- 
presentation of all possible trajectories starting at a collision point within a 
small domain around this point. 

4 

Subsequently, Sundman (1907-1912) made the second and decisive con- 
tribution to the analytical solution of the general three -body problem. He 
demonstrated the possibility of representing the coordinates and the veloci- 
ties of the three bodies by means of convergent series and determined their 
radius of convergence. The key points of his proof are: a) the existence 
theorem by Cauchy- Picard for the solution of a system of differential equa- 
tions of the 1st order; and b) the transformation of the independent variable 
t into a new variable u by means of 

(I. 1) du = “ dt 

where r is the distance between two colliding bodies, that is the distance 
which will vanish at a given value t of t. After showing that lirn u exists 
and is finite, Sundman demonstrated that the coordinates and vetocities, 
even in the neighborhood of a binary collision point t= t, can be expanded in 
convergent power series of (t - t) 1 ' 3 . This helps the understanding of the 
analytical nature of the singularity point. In fact, three branches of the 
same function can be permuted one into another around the point t = T. It 
is possible, therefore, to continue analytically the representation of u, the 
coordinates and the velocities after, or before, the collision point by always 
taking the real value of the cubic root of t - t . 
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It is of historical interest, however, to note that in order to remove the 
singularity at a binary collision point t = T, the appropriate choice of the 
variable u, according to th^ transformation (I. 1) had been pointed out earlier, 
but not exploited, by Bruns (1884). , 

After Sundman's work, the regularizing variable concept was generalized 
in various directions (successive binary collisions, triple collisions, imagin- 
ary collisions, other regularizing transformations, interpretation in the com- 
plex domain, etc. ) and also gained acceptance among authors of classical 
books on Celestial Mechanics, notably, Charlier , Wintner?, Siegel® and 
others. Recently, Arenstorff^ has published a series of papers on this sub- 
ject using methods taken from the theory of the functions of complex variables. 
The goal of this paper is to describe an algorithm leading to the construction 
of convergent series expansions of the coordinates when a regularizing vari- 
able of Sundman's type is used as the independent variable. Our choice for 
this variable is given by the transformation 

(I. 2) du = Vdt 

indicated by Levi-Civita at the end of his celebrated memoir of the year 
1917. ^ In (I. 2) V is the total potential function defined as the negative of 
the total potential energy. Vernic^ advocated also the application of the 
transformation (1.2) to numerical computations. A goal similar to ours is 
that pursued by Steffensen^ to obtain the time series expansion in the case 
of the planar restricted three -body problem by means of recursion formulas. 

It turns out that replacing t by u, recursion formulas can also be found for 
the case of the general three-body problem. In fact, explicit expressions for 
the derivative of the coordinates dx. , for sufficiently large values of v will 

IhF 

be established later in this report. The application of the transformation 
(I. 1) or (1.2) to the 2-body problem leads to a unified set of formulas valid 
for any kind of conic sections. A heuristic approach to the derivation of such 
universal formulas has been illustrated elsewhere. ^ 


II. THE EQUATIONS OF MOTION AND 
THEIR SOLUTION BY TIME POWER SERIES 

Let niQ, m ^ , m^ be three non-zero masses and write the equations of 
motion of all three bodies in an inertial Cartesian reference frame 


(II. 1) 


av 



(i = 0,1,2) 
(x "*y “* z ) 
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Here, the dots represent derivatives with respect to the variable t 

(II. 2) T = kt, 

k = gravitational constant, 
t = time. 

The total potential function V in (II. 1) is defined by 

(IX. 3) V=Y m.m. — , (i * j ) 

. . i j r. . 

1J J 

where the mutual distances r.. between bodies are given by 

1J I 

r 2 2 2 ~|2 

(II. 4) r„ = [Jx. - Xj) + <Yj - y.) + (Zj - z.) J 

Then, the explicit expression for equation (II. 1) is 
(II. 5) x. = m. ^ . (x. - x.) + m k M - x. ) . 

where the subscripts i, j, k are permuted cyclically according to the 
following rule 


i 

j 

k 

0 

1 

2 

1 

2 

0 

2 

0 

1 


and, in general, 


a. is the inverse cube of the distance r.. 


(II. 7) 


1 



Equations similar to (II. 5) can be written for the other coordinates y. and z^. 


The formal power series solution of the equations of type (II. 5) can be 
written as follows 

» (i = 0,1,2) 

(II. 8) *.co W ' 

V =0 


(x. -* y. — z.) 

' l 7 i l 
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where 

(II. 9) 


d V x. 


X iv v! [ v] t = 0 
d t 


and t- 0 is an initial (origin) value of T. The first two coefficients of the 
expansion (II. 8) are 


(II. 10) x = x. (0) , x = x. (0) , 

i0 i ll l 


that is, position and velocity components at origin t= 0. Initial position and 
velocities components of all three bodies are assumed to be known. 

The next coefficients x. ^ of the Maclaurin expansion (II. 8) is — x. (0) 
which can also be obtained 1 by directly evaluating the RHS* of equation (II. 5) 
at T = 0. The other coefficients of the higher order terms can be obtained by 
repeated differentiation with respect to T of equation (II. 5). 

This repeated differentiation requires, however, a considerable effort 
in algebraic manipulations. In fact, there are nine differential equations of 
2nd order, coupled together by the mutual distances; hence, there are nine 
Maclaurin series expansions in which the coefficients depend on several in- 
termixed parameters (mutual distances, masses, positions and velocities). 
The expressions for these coefficients become extremely complicated, and 
their complexity increases with the order of the differentiation. This ex- 
plains why in the past the computation of these coefficients was considered 
to be impractical, until Steffensen found a relatively easy way of computing 
them in a special case. This author showed that it is possible in the planar 
restricted case of the three-body problem, which reduces to only two coupled 
2nd order differential equations, to establish recursion formulas for the 
computation of the coefficients of the two corresponding time series. Earlier 
attempts on this sublet were those for the asteroidal case of the three -body 
problem by Stumpff and Sconzo ^ . Rabe^, Deprit^ et al. have recently 
made use of Steffensen' s formulation in numerical studies of planar restricted 
problems. Another interesting attempt to get the high order terms of the 
time series has been described by Grobner*® who applied Lie's concept of 
infinitesimal transformation^ to the study of the n-body problem. In the 
application of this concept to the planar case of the three -body system formed 
by the Earth, Moon and spaceship, however, only the 4th order coefficient of 
the time series has been computed explicitly. In another paper by 
Bogayevskiy^^ , where a lengthy recursion formula for Grobner coefficients 
is presented, only the 3rd order coefficient is given explicitly in terms of the 
initial conditions. Bogayevskiy 's procedure is complicated not only because 
the notation used is very complex but also because the number of terms which 
*Right-hand side 
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constitute the desired expression is to be determined solving Diophantine 
equations. 


Now, we deem it noteworthy to present the expressions for x and x q4 
arrived at by two successive differentiations of equation (II. 5) written for 
i = 0. The following concise expressions can, in fact, be derived 


(II. 12) 


3 

< n - U > X 03 = l! [ ] = T[ m l W 0i a 01 +m 2%2 a 02] (o) V O) 


+ m 2 M 02 J 


( 0 ) 


( 0 ) 

V°> 


‘ 6 [ m l M 01 

- 1 [ m i w oi a oi] (0 x i <0) + i[ m i t ‘oi] 0 x i (0) 


-K m 2%2 a 02] (0) X 2 (0) 


+ i[ m i^ 


02j 


( 0 ) 


x 2 (0) , 


d 4 x 

= Tl [ ] (0) = [ 8( m l M 01 ( %l' 5CT 01 ) + m 2 M 02 (e 02- 5a 02 ) } 


04 


dT 

+ Ti { (n V 0l +n, 2** 02 )Z + m 0 (m 1% 1 + m 2 M 02>}] X 0 (0) 


+ i[ m l /J 0i a 01 + m 2*‘02 a 02] (o X 0 (0) 

+ m l [- boi <*01 - Hn> + l4{ m 2^02^12- (m 0 +m l );i 0 2 l 

- m 2^ 01 ( %2 + ^2>}] (0) X l (0) - i[ m i^ 0 i a 0l] (()) i l <0) 

+ m 2 [-¥%2 (t 02- 5a 02> + TT{ m i%l M 12 - (m 0 +m 2^02 

' m l M 02 (M 01 +Ai 12 ) }] (0 X 2 (0) ‘ 4[ m 2 A< 02 a 02] X 2 <0) * 


if the following notations are adopted 
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(11.13) 

(11.14) 
where 


a = — s 
01 2 01 , 

r oi 

€ oi = T s oi ’ 

r 01 


s oi = 


m. 


01 


(y S z) [ (i r*0 )2 ‘ - m 2 <X r X 0 ) { M ffi (x 2' X 0 ) - M 12 (x 2- X l ) }l 


01 


The symbol S X. . . J means a sum extended to terms in the variables 


(y.z ) 1 

y and z similar to that in x. 


Inspecting the structure of the RHS of both equations (II. 11) and (II. 12) 
we see that up to the third order term the perturbation effect upon the co- 
ordinates of the three bodies is of the first order with respect to the masses, 
while starting from the 4th order term this effect becomes of 2nd order or 
higher. 


Applying the cyclic permutation of the indices indicated by (II. 6) to 
the equations (II. 13) and (II. 14), the expressions for x , x and x , x 
can easily be derived from (II. 11) and (II. 12), respectively. Similar 24 
expressions for the y and z coordinates can also be easily derived. We ob- 
serve that Bogayevskiy's expression for x , presented in ten strings of 
terms, can be simplified to only 8 non-vanishing terms and in its simplified 
version it coincides with our polynomial expression (II. 11) which contains 
precisely 8 terms. 


We want now to call the attention of the reader to the fact that both equa- 
tions (II. 11) and (II. 12) have been brought to a form analogous to that of the 
coefficients of the Lagrangian f and g series in the two body problem. The 
symbols n , a.., c .. are also an extension of the symbols ju, a, « adopted 
for the two -body problem in the referenced paper. Although the said analogy 
can be extended further and formulas similar to Cipolletti's^ recursion for- 
mulas could be established (this will be the object of a later IBM funded in- 
vestigation) here, we will not try to derive explicitly the expressions for 
the derivatives of the coordinates with respect to the variable t corres- 
ponding to an order v ^ 5. Our aim as it has been stated in the introduction 
is to obtain these high order derivatives with respect to the new variable u. 
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It is important, however, to determine the radius of convergence of the 
series expansion (II. 8). We dedicate the next section to this problem. 

HI, CONVERGENCE OF THE TIME SERIES EXPANSION 

It can be said, in general, that the convergence radius of the series 
expansion (II. 8) is determined by the distance of the nearest singularity 
point (real or imaginary) to the origin. It can be demonstrated that the 
series (II. 8) is convergent in the neighborhood of a non- collision point. An 
upper positive bound T of the independent variable T reckoned from the 
time of a non-collision point T can, in fact, be found such that within the 
interval (-T, T), centered at T, the series (II. 8) is convergent. To 
demonstrate this we will use Cauchy's theorem as implemented by Picard^ 
on the existence and uniqueness of the solution of a system of first order 
differential equations in the real domain. 

Let 

(III. 1) f i [ x 1 (T)> x 2 (T). x n (T) ]’ (i = 1 . 2, .... n ) 

be real functions of the real independent variable T , but not containing T 
explicitly, satisfying the following conditions 


a) they can be expanded in convergent power series of x^- x^ 
when 

(HI. 2) | x. - x. | < Tl. . 

where x. = x. (T ) , T is an arbitrary value of T, which can be taken T= 0 F and 
all T|^ are positive quantities; 

b) it is 

(III. 3) I L I < Fi . 

when x. satisfies the inequality (HI. 2) and where all are positive quantities. 
Then, ihe system of differential equations 


(in. 4) 



admits one and only one solution such that 
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(III. 5) lim x. = x. . 

T-r i i 

Under the conditions specified above it can be added that it is possible to 
expand the unknown functions x. (T ) into a convergent power series of T - T* 
for any value of T such that 

(III. 6) I T- T I < T , 

where 

(HI. 7) 

Furthermore, the inequality (IIL 2) is verified when T is chosen according to 

( 111 . 6 ). 



In order to apply the Cauchy-Picard theorem stated above to our system 
of nine differential equations of 2nd order (II. 5), we transform this system 
into another system of 18 first order equations as follows 


(in. 8) 

dx. 

i 



dT 

X i * 


dx. 

<W z i ) 

(HI. 9) 

i 

dT 

= ™/ u (yV + I V i ik (x k- x i ) - 


Now, let T|q and T] be two positive constants such that 

(in. 1°) |x i -x.|, ly.-yj, 

1 v 

A 

o 



U = 0,1,2) 

(in. 11) |x.-x.|, |y.-*.|. 

K - 

A 

o 


and we will demonstrate that condition a) of Cauchy's theorem is satisfied, 
that is the RHS of equations (III, 8) and (III. 9) can be expanded in convergent 
power series of the differences x. - x. , etc. It suffices to demonstrate that 
under the conditions (III. 10) the inverse cube of any of the three distances 
r.. can be expanded into a convergent series. 


We take into consideration, for example, the distance r and we write 


(III. 12) 


01 


-2 
= r oi 


+ P 


' 01 


where 


105 


TRAJECTORY ANALYSIS AND GUIDANCE THEORY 


(III. 13) 


_2 2 /JFi 

r oi = r oi (T) 


and P is a second degree polynomial in the differences x Q - x Qf etc. 

The explicit expression of ^ is 

p oi = (f.z) [ 2( *rV{ (x i‘ 3E i ) ‘ (x o' s o ) } + { (x r s i ) ‘ (x o' 5 o ) } ]• 


It follows that 
(III. 15) 
because 

(III. 16) 


l^l 1 < 12 Vo + 12l ^o 


|x r x 0 l , ly x -y 0 l. L,- 0 I « r oi • 

Consequently, the expansion of — as well as that of W ^ 

r oi 


(III. 17) 


01 


1 , 01 v 
— (1+ — > 


3 

2 


01 


01 


01 


is convergent if we choose T] such that 


(III. 18) 


l2 *oi't + T oi 


(III. 19) 


^0 < 


01 


6+ 4/3 

Condition (III. 19) is satisfied a fortiori if we take 
(III. 20) 


01 


^0~ 14 * 


We determine now two upper bounds for the RHS of both equations 
(III. 8) and (III. 9). We begin with the equations of type (III. 9) and we observe 
that with the choice (III. 20) we can first deduce 


(III. 21) 


r oi > / 


4 • i2F oi% - u ”- 


7 r 01 
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and since 


Cl - x 0 U I x ! - x 0 I + l x o - V + I X 1 - x il . 


we deduce also 
(III. 22) 


I X 1 - *ol < ? 01 + 2ll o = f f or 


Then 


(m.23| Vl*(-Sr/f 


8 _ 

r 


01 , 2 
411 0 


and other inequalities similar to (III. 23) can be established for 

•%2 (X 2 - X 0 ) l* etc - 


Now, if T - T is a non-collision point we can find a lower bound of the 
three distances at time T = T , that is a positive number r}^r)QSuch that 


(IK. 24) 


r 01’ r 02’ r 1 2 


> 


1 4r| . 


Thus, the inequalities (III. 10) can be replaced by 


<m.25) |x. - x.|, |y. -yj, | z. - z. | < t, 

and the inequality (III. 23) by 


(III. 26) 


*oi (x i - V* < 72 

4*n 


Considering the equation (III. 9) we then obtain 


(U1.27, I 3^1 < 7 < —z ■ 

4r] J 4r) 

where 

(III. 28) M = m^ + nij + . 

Now, we deal with the equations of type (III. 8). We first deduce the 
existence of upper bounds for the velocities at time T= T. We prove this 
by considering the energy integral 
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(III. 29) 


1 V f Z X * 2 -L - Z \ 

2 1 (x i + y i + Z i > - 

i= 0 


V = h 


rewritten as follows 
(III. 30) 

where 

(III. 31) 


i= 0 


m. (x 2 + y? + z 2 ) - 2 

- l l i 


m 0 m i nl 2 „ m 0 m i m 2 u 

U = h , 


M 


M 


m 0 m l m 2 _ _ m Q m i + m 0 m 2 + m i m 2 


M 


01 


02 


12 


and h is a constant. 

When T-»? , U -* U, where U according to (III. 31) and (III. 24) satisfies 
the inequality 


(IU. 32) 


m 0 m i m 2 - 


M 


U < + na Q m 2 + m^) 


Now, we can apply to the RHS of (III. 32) the following obvious result 
M 2 = (m Q + + t" 2 ) 2 > 3(1^0^ + m Q m 2 + 

Doing so, the inequality (III. 32) becomes 


(III. 33) 


rri m , m_ . 

0 1 2 jj 4 M 


M 


42 T) * 


It can also be proved that, 


M > 4m Q m 1 , 4m Q m 2 , 


then, we have in general 
(III. 34) 


w m.m 

M l k 


it > -L 

4 M 


(j 4 k) 


If we now divide by m all terms of equation (IU. 30) written for T = T, where 

m = min [ m 0 , m^ , m^) , 
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we deduce 

Letting 

(in. 35) 

we have 


t2 t2 tZ 

x. , y. , z. < 


M . M. ... 

21mr| 4 * h ' 


f, = /2T 

'o V21mri 


+ -tM 


z. < T| _ 

l 1 '0 


1^1 • I yj. 

and we obtain from (III. II) 

(III. 36) |x.| f | y. |. |z.| < 2 r| 0 . 

We may conclude that selecting T as follows 


(III. 37) 


T = min 


r n. o \ 

l2f| 0 ’(M/4Tl 2 ) J * 


the series expansion (II. 8) is convergent for any T such that 
(in. 38) | T - T | < T. 


The entire proof given above has been borrowed, with only slight modi- 
fications, from the original memoir by Sundman 4c . It is worth noting that, 
of both ratios contained in the RHS of equation (in. 37), the first is smaller 
than the second. In fact 

- , 2 

Tl , ^0 T] ( 8T1Ti 0 , \ r\ ( 8M , \ 

■ 2 ^ (mTSti 2 ) = 2 ?^ ( — ' = W Q +2T lN-l)> o 

because 3m ^ M, and consequently > 1, > 1. 

7 24m 21m 

Hence, all the series expansions of type (II. 8) are convergent for any T 
such that 


(IU. 39) 


|T- T | < 


l 


4M 
21m T| 


+ M | h | 
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and these series represent the analytical solution of the three-body problem 
which satisfies the preset initial conditions for position and velocity at time 
T «T. 


IV. ON REGULARIZING VARIABLES 


We say that a variable is a regularizing variable if it has the property 
of removing that singularity of the differential equations of the motion which 
corresponds to a binary collision. 


We will show that the variable u introduced by the Levi-Civita trans- 
formation (1.2) is a regularizing variable. We demonstrate first the fol- 
lowing lemma: if S(r . r , r ) is a symmetric and homogeneous function 
of first degree of the three distances r Qj» r 02* r 12* ^ ien var ^ a ^^ e u 
defined by means of the differential operator 


(IV. 1) 


d _ 1_ d_ 

dT S du 


is a regularizing variable. 

In fact, applying twice the operator (IV. 1 ) to x. and denoting by primes 


the derivatives with respect to u, 

we obtain 

(IV. 2) L = 5 x i' 

„ 1 ** 

S' , 

(IV. 3) x. = - 2 x. - 

’ 

or vice versa 

(IV. 4) x' = s X. 

» 

. 

0 2.. 

(XV. 5) X~ = SSx. + 

S X i ' 

if we observe that SS = S'. 

By virtue of (IV. 3) and (IV. 5), the equation of motion (II. 1) acquires 

the form 

/> s' * , 

sf. 8_V_ 

(IV. 6) X. = gX.+ 

m. 9 x. 

or the equivalent form 

l i 
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x."= SSx. + £ 

i i m. 9x. 

x i 


(IV. 7) 

The latter can also be rewritten as follows 


(iv. 8) x ~ * si. V 20L V JL5L ri J 

x x Z_..9r ij m. Z^.Or.. 9x. 

i ij xj x 

We examine now the order of the infinitesimal terms which constitute 
the RHS of equation (IV. 8) written for i = 0, 1 when the two corresponding 
bodies collide at a certain time T = T . that is when ^lim r Q1 = 0. Considering 

r 01 as an infini t e simal quantity of first order, the order of all functions de- 
pending on r can be listed as follows 


Function 

Order 


x. (and also 
1 \ ) 9r 


as 


01 


01 


0 


f 2 


av 

ar oi 

-2 


9r 


01 


a V a Xj 


All these results are evident except those for x. and r which may be de- 
duced as a consequence of the following four limits known as Bisconcini- 
Sundman's^’ relationships 


U -7 (/^ f 01 ) = H = - / 2 < m o +m l ) 


(IV. 9) 


lim/ 


lim 
T -*T ' 

lim , 


where C^, C^, are three constants such that 

c 2 + c 2 + C 2 = 1 

X y z 


k.) 

|= C H 

l) 

X 

Yi 

) = C H 

X 

' Y 

z. 

)= C H 

X 

/ z 


(i = 0,1) 


Then, it follows from (IV. 8) 
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where 


lim 

a = - u 

T T 


We conclude that 

(IV. 10) lm L X* = c it 0, (i, = 0,1) , 

\ 1 ~* U 1 

that is, the RHS of equation (IV. 8) is a constant C different from zero when 
the two bodies (i = 0,1) collide. 

Asa consequence of this lemma we conclude that the variable u defined 
by Levi-Civita differential operator 


(IV. 11) 



is also a regularizing variable. 


V. USE OF LEVI CIVITA'S DIFFERENTIAL 
OPERATOR AND RECURSION FORMULAS FOR THE SOLUTION 


Applying the differential operator (IV. 11) to the equation of motion 
(II. 1 ), we obtain 


(V. 1) 


Vx. + Vx. 
l i 


i av 

m.V 8x . 


(i = 0,1,2) 
<V y i~ z i> 


In this equation the function V is to be considered as a function of the variable 
u. Since V is a known function of T, in order to obtain its dependence on u 
we proceed as follows. First, we calculate 

T 

(V.2) u ( T ) = \ V(5) d§ 

J 0 


The RHS of equation (V. 2) is a power series in T , which for practical pur- 
poses can be truncated to a polynomial of certain degree. 

Now, reversing the series given by (V.2), we get 
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(V. 3) t = t (u) ; 

* 

then we insert (V. 3) into the expression of V(T). Denoting by V(u) the 
result of this operation we can rewrite equation (V. 1) in the new form 

(V. 4) V x . + V x = X. , 


* i *-i dv 

X. = — V . 

i m . 8x . 

l l 


*i (u) = I *iv u " 


be the formal series expansion of the solution of equation (V. 4). Letting 
also 


* r* * v 

X i<u) = I X iv u * 


* * V 

V=y V v a 


and inserting (V. 6), (V. 7) and (V. 8) into equation (V. 4) it will not be difficult 
to arrive at the following recursion formula 


1 v - 1 

'-1)^. X iV-2 vv. I ^ ' k) V k x iv -k' < v ^ 2) * 


% j)c 

Formulas similar to (V. 9) can also be established for y. z.^ after defining 

(V.10) f = J_ V • 1 -|i , 

l m. 9y. 


* 1 *-l 8V 

Z. = —V 4 

i m . 9 z . 
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(V. 12) 


x i0 = X i ( °> 

(V. 13) 


* 1 
x il " -5- *i<°> 
v o 

are known, the use of (V. 9) for v = Z, 3, . . . provides all the coefficients 
needed to write the formal expansion ^f x. in a^power series of u. We need 
to start from known values of X. Q , V Q and V as functions of the initial 

conditions. For this purpose we observe that starting from the given initial 
conditions (II. 10) we can write the following linear approximation 

(i = 0,1,2) 

(V. 14) x. = x. + x.T, . 

' ' i lO il (Xi -y. z.) 

Then, we have 


(V. 15) 


2 _2 T 
r ij " r ij + P ijl * 

where 


7 ij = (y s )Z )[ x jo- x io] • 



P ijl= 2 (y, z) [ (X jO ' X iO^ X jl ' X il * ] ' 

It follows from (V. 15) that 
m.m. 

(V . 16) _i_,L = + g T . 

where 


m.m. 

g ij0 " 

ij 

Hence 


g ijl ” ’ 8i i° P ijl ' 

ij 

(V. 17) 


v = v 0 + V, x , 
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where 


0 = I. g ijO ’ V 1 = I. g ijl * 
l J U 


and according to (V. 2) we obtain 
(V. 18) 


U(T ) = V 0 T + j VjT 2 . 


Reversing (V. 18) we get 
(V. 19) 


V 

1 112 
= T U ■ 2 “3 u • 
V 0 


Inserting (V. 19) into (V. 17) and truncating the result to the linear approxi~ 
mation we obtain 

# * # 

V(u) = V Q + VjU , 

* 

*-l 1 V l 

v <«)- — - *£ « . 

v o v o 


where 
(V. 20) 


* 

V„ 


v 0- 


* V 1 
V = — - 

1 v o ' 


We have also 
(V. 21 ) 


1 

"ij = IT 

r. . 
1J 


2r. 5 P « 1 

ij 


(V. 22) 


x j ■ x i = (x j0 - x i0> + (x ji - X il> T 


Finally, inserting (V. 19) into (V. 21) and (V. 22) and performing the products 
“ x.)V , x^)V and their sum we obtain 

* * * 
x i = x i0 + x n“ • 
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where X. and X are now known expressions of the initial conditions (II. 10). 
As a matPer of i:act we can find directly from its definition 


* 

X i2 



Applying twice the inverse operator of (IV. 11) to x i we have in fact 


1 d S 11 d/1^.1/1 !l!Li . I dV ^J_\ 

2 = 2 V dFV V dr ) - 2 V y 2 dT 2 y 3 dr dr ) 


d 2 X! j dy dx. 


and the evaluation of both sides at the origin (u = 0 and T = 0, respectively) 
provides 

1 

(V. 23) 


£ = — — [? ] 
i2 _ ,.2 L Vt=0 


2 


2 V 0 


3 V l X il 


where 




t=0 


The result expressed by (V. 23) coincides with that obtained from (V. 9) 
if we put V= 2. We take into consideration (V. 13) and (V. 20), and we ob- 
serve that by virtue of (V. 5) * i0 = v KiW 

£ o £ * 

Next, in order to obtain x , we need to know V 2> We arrive at V 2 by 
replacing (V. 1 4) by a second dfegree polynomial in u which is known because 
x. is now a known quantity. The computation of the 2nd order expressions 
in 2 the variable u for r.., m.m. and V will provide the desired The same 

r . . 

procedure is then applied to compute V^, .... and, therefore, all other co- 
efficients of the series expansion (V.6). 

The algebra involved in this computation, in summary, consists of the 
following operations to be performed on polynomials 

. subtraction 
. multiplication ^ 

. integer and - ^ and - — power 
. reversion and inversion 
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• substitution of the variable of a polynomial by another polynomial 
. integration 

This algebra is being programmed for polynomials of arbitrarily high 
degree. 

We observe that the polynomials in the variable T (polynomials of first 
degree) of type (V. 14) are needed only to initialize the computation. The 
computation is then continued by using polynomials in the variable u. At 
each successive approximation, a new term will be added to the polynomials 
expressing the nine coordinates. 

The relationship between T and u should also be computed at each ap- 
proximation replacing equation (V. 2) by the equivalent equation 

r u 1 

(V. 24) T = [ da . 

0 V(a) 


CONCLUSION 

The series expansion of the solution of the three-body problem has been 
found in terms of Levi- Civita's regularizing variable u. The recursion for- 
mulas of type (V, 9) are the key points of the computation. By means of these 
formulas the coefficients of the series expansions for all nine coordinates 
are expressed as functions of the masses and the initial conditions. 

That these series are convergent can be demonstrated using a well 
known theorem applied to the solution of linear differential equations of type 
(V. 1) when they are solved by the power series method. This proof is being 
adapted to our type of equations and it will be presented in the next report. 

It is expected that by using power series in the variable u, the analytical 
representation of the motion can be extended to an interval of time much 
larger than the step- size used in numerical integration procedures. If so, 
the method could advantageously be utilized in place of numerical integration. 
It would, in fact, lend itself to "being exploited as a numerical tool for the 
analytical continuation of the solution. A previous effort in this direction, 
using time series expansions, was made in 1955 by the author of this inves- 
tigation. 

In the next report we also intend to present a computed numerical case 
of the three -body problem using the method illustrated above. Specifically, 
we will choose an example considered as a "zwackmassig 11 one by Bohlin ^ 
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and for which results, both in tabular and graphical form, have been obtain* 
by Zumkley who used a numerical integration procedure. In this examp 
the masses and the distances have the same order of magnitude, and the 
example is, therefore, appropriate to be used as a test case of the method 
described above. The only simplification used in this example is that it 
deals with a planar case instead of a tridimensional one, and this has been 
done in order to reduce the burden of the computation without prejudicing th 
general validity of the formulation. Some preliminary computed results we 
have obtained by hand computation, carried out up to 5th order terms in u, 
confirm the expectation that large intervals of time can be covered by this 
kind of truncated series representation of the motion of all three bodies. 
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ABSTRACT 

By utilizing results of Hamiltonian theory and the von Zeipel 
method for treating artificial satellite orbits, error bounds are derived 
for a general class of orbits with eccentricity less than one. In order 
to extend the error bounds for the general axisymmetric problem to time 
intervals of the order l/Jg, the known integral of energy is utilized 
to calibrate the governing differential equations for the rapidly rotating 
phase. The non-singular rapid phase in this analysis is taken to be the 
sum of the mean anomaly, argument of periapsis and the right ascension of 
the ascending node. A corresponding analysis for the general asymmetric 
problem (including the tesseral harmonics) is also given. From the 
general error analysis an algorithm is derived for the computation of 
the correct initial conditions consistent with the expected accuracy 
of the theory. Numerical results verifying the conclusions of the theory 
presented in this paper are also given. 

' This work was performed in association with research sponsored by the National Aeronautics and 
Space Administration under Research Grant NsG-133-61 
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I. INTRODUCTION 

The analytical theory of artificial satellite motion has been the 
subject of very intensive study since the launching of the first 
artificial satellite in 1957. In fact, many aspects of the problem 
had been studied before that time in connection with the theories of 
celestial mechanics. The result of the study has been a very extensive 
list of papers offering solutions of many differing forms and techniques 
of achieving them. However, with the exception of the work of Kyner 
(Ref. 1), no other solution is known to the authors that offers rigorous 
error bounds on the position and velocity for a general class of orbits, 
e.g., inclined orbits of any eccentricity less than one. Naturally, the 
orbits at critical inclination and orbits in resonance with the tesseral 
harmonics must be excepted from the general class. It is then a matter 
of general interest to derive such error bounds. 

From a fundamental point of view, the problem of artificial satel- 
lite motion can be classified as a special case of a general class of 
non-linear oscillation problems. Non-linear oscillation problems can 
be treated with varying degrees of success by the general averaging 
methods developed by Krylov, Bogoliubov and Mitropolskii (i.e., Ref. 2). 
For these methods of averaging there exists an associated technique for 
establishing bounds on the error build-up in a specified time between the 
exact and the approximate solutions (first order or higher order). Now, 
the method of application of the technique of averaging to the problem 
of artificial satellite motion depends rather heavily on the particular 
choice of variables employed. In the case of Kyner’ s work, averaging 
could be applied directly; in most other approaches to the problem the 
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use of averaging is more or less disguised. 

One of the most widely used perturbation methods in treating 

artificial satellite orbits has been the method of von Zeipel as adopted 

by Brouwer (Ref. 3) and Kozai (Ref. 4). This method is one of successive 

canonical transformations and is necessarily carried out in the variables 

of Delaunay (L,G,H,^,g,h) . With a slight change of variables and a 

choice of a different intermediary orbit, the same method was applied 

by Garfinkel (Refs. 5,6). Furthermore, it has been shown (Refs. 7,8) 

that the von Zeipel method of canonical transformations is a particular 

form of the method of averaging. Hence, by drawing on the equivalence 

to averaging, rigorous error bounds could be established for the 

Delaunay variables directly. Unfortunately, bounds obtainable in this 

way for the Delaunay variables £ and g are unsatisfactory for very 

small eccentricity (i.e., e' < J where J 0 is the oblateness 

z 2 

—3 

parameter of 0(10 )) due to a singularity at zero eccentricity in the 

short period terms. A further drawback is the singularity at zero 
inclination. Since no singularities exist in the coordinates for zero 
eccentricity and/or inclination, one would expect that these objections 
to the bounds would not exist for a suitable choice of variables. The 
error bounds derived by such direct application of differential equation 
theory turn out to be unsatisfactory for large time intervals i.e., time 
intervals of the order l/J^ " Since one of the problems of interest in 
applying closed— form orbit theories is orbit prediction over long periods 
of time, the error theory must be modified. The modification is a more 
involved problem and a separate treatment is presented here. 

In this report, the problem is analyzed in canonical variables; 
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the three sets of interest are those due to Delaunay, Hill and PoincareC 
Of these variables, the Poincare^ set is non-singular for both zero 
eccentricity and inclination, the Hill set singular for zero inclination 
and the Delaunay set singular for both zero inclination and eccentricity. 
The advantages of the Hill set are the simple forms of the in-plane 
coordinate perturbations which are obtained directly from known 
generating functions. It was shown by Izsak (Ref. 9) that, to first 
order in the oblateness coefficient , the in-plane position and 
velocity components of a satellite are obtainable by converting via 
Keplerian formulae from Brouwer's averaged Delaunay variables (L 1 ,G' ,H' , 
j}',g') to corresponding "averaged” position and velocity and then 
superimposing the short— period fluctuations. These short— period fluctu- 
ations were shown to be obtainable by rewriting Brouwer's short-period 
generating function in terms of the Hill variables and taking 

appropriate partial derivatives. These short-period fluctuations are 
well-behaved (unlike those in i,g) when eccentricity goes to zero. 

Recent investigations by Vagners (Ref. 10) have obtained in the same 
manner first order long-period fluctuations in the Hill variables by 
rewriting Brouwer's long-period generating function relating (L ( G , 

H” , i” , g" , h”) to (L* ,G' , H ' t V , g* ,h* ) , including general formulas for 
the effects of any zonal harmonic. Analogous "medium-period" (i.e., 
daily) fluctuations in the Hill variables were obtained in a general 
form for the effects of the tesseral and sectorial harmonics. Since the 
analysis given by Vagners was applicable to any set of canonical variables, 
then similar results could readily be obtained for the Poincare'" variables . 

Utilizing the results of Izsak and Vagners, an analysis is carried 
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out in this report which parallels every canonical tranf ormation of 
the Delaunay variables by an appropriate canonical transformation of 
• some general set of canonical variables including the removal of second 
order short-period terms from the Hamiltonian, In this way, rigorous 
error bounds on the first-order solution are established which are 
independent of the eccentricity for Hill variables and independent of 
eccentricity and inclination for the Poincare / variables (as long as e 
is not too close to one). As is shown, these bounds are unsatisfactory 
for long time intervals and another method is offered, 

A discussion is presented of the various terms arising in the error 
bound. Particular attention is focused on the question of initial 
condition errors; this question is of interest when computing by means 
of a "closed-form" satellite theory a satellite* s ephemeris from some 
given initial position and velocity vectors. In view of the extensive 
comparison studies of different orbit theories conducted by Arsenault, 
Enright and Purcell (Ref. 11), wherein the problem of initialization 
plays such an important role, this question assumes considerable 
importance. An energy method is then given for greatly decreasing the 
primary in- track position error build-up due to initial conditions and 
some typical results are quoted. The algorithm of computing the correct 
initial conditions arises directly from the extended error bound theory. 

The authors wish to acknowledge the contribution of Small (Ref. 12), 
who first utilized the energy method in reducing initialization errors 
in his solution to the problem of satellite motion about an oblate planet. 
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II. GENERAL BOUNDS ON SATELLITE MOTION 


Before proceeding to more specific treatment of the error problem, 
some general statements concerning the a priori bounds on the motion may 
be made. First, one can consider the motion of a satellite in a general 
axi-symmetric gravitational field for which two integrals of the motion 
are known. If the potential field is represented by 



oo . 

\ N 

f- 

1 

II 

> 

1 - I J N (■ 

N=2 

r) V sin P> 


[ U (r) + U x (r ,p)] (II-1) 


where [3 is the latitude, the equatorial radius, p. the gravita- 

tional constant, r the radius and J N numerical coefficients, then 
it can readily be shown that the total energy and the polar component 
of the angular momentum are constants of the motion. The two exact 
integrals may be written in the form 


i + i =h 


and 


= v/pad - e 2 ) cos i = k 


(II-2) 

(II-3) 


where a is the semi-major axis of the orbit, e the eccentricity, i 
the orbital inclination and k^ k 2 are constants. 

The two integrals (II-2) and (II-3) imply that if k 1 and k 2 
are given, then the motion of the satellite is confined to a region 
bounded by a "zero velocity" surface (Ref. 13). With initial conditions 
specifying and k 2 one can write the a priori bounds in the form 


0 < &i( k i l k 2 , 6 ) <r< 6 2 (k l ,k 2 ,G) 


(II-4) 


where e = J 2 , k 1 > 0, k 2 > 0 and are assumed values of 0(1). 
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General bounds of this type are developed by Poritsky (Ref. 14) and 
given for € = 0 by Kyner (Ref. 1). Here the explicit forms of & 
and Sg are not of direct interest. 

In the more general problem of a longitude dependent potential one 
no longer has the two integrals (II-2) and (II-3). Such a potential 
arises when one includes the tesseral harmonics of the Earth's field in 
the general satellite problem. However, by considering a rotating coor- 
dinate system fixed in the primary, one can readily determine the Jacobi 
integral of the system. In this case one specifies only the upper bound 
by the zero- velocity surface. 

One assumes then that a priori bounds on the state vector x are 
known; namely, if the initial state vector x(0) is in a set D, then 

|x| < C(x(0) , e) (II-5) 

where the solution depends on a small parameter e . Since for near- 
earth satellites one is concerned with elliptical orbits, the set D 
will be specified by the requirement of negative energy and a non-zero 
initial value of the angular momentum. If the state vector chosen for 
the description of the motion is some canonical set (q,p), then the 
equations of motion take the form 
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Then, since ^ ~ is continuous and satisfies a Lipschitz condition 
locally in x in some bounded region (then) a solution for all t 
exists as a consequence of (II-5) . 

Note that implicit in (II-5) is also a restriction on how close the 
energy and the angular momentum may be to zero. For the general bounds 
to hold, these initial values must be sufficiently different from zero 
so that the perturbations, of order e in the satellite problem, do 
not cause the state vector x to become arbitrarily large. 
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III, THE SECOND-ORDER HAMILTONIAN 

Inherent in a specific discussion of error bounds is a knowledge of 
the characteristics of the analytical method used in the fundamental 
solution and a knowledge of the behavior of various functions arising 
therein. The method utilized in the following analysis is the von Zeipel 
method and the system analyzed is a Hamiltonian system. For a brief 
review of the von Zeipel procedure, the reader is referred to Ref. 10; 
the specific details of the orbit problem solution may be found in 
Refs. 3 and 4. 

It turns out to be convenient to introduce the three sets of 
canonical variables due to Delaunay, Hill and Poincare' . (Recall that 
the original solution of Brouwer was carried out in Delaunay variables.) 
These sets of variables are defined in the following manner: The 

Delaunay variables, denoted by y, are given as 



where £ is the mean anomaly 

g = w the argument of pericenter 
h = H the right ascension of the ascending node 
L ss yj [_ia 

G = y / (1 - e2) 

H = G cos i 
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Denoting the Hill variables by z, one finds 


r 
u 
h 
R 
G 
_ H 

where u the central angle or argument of latitude 
R = r the radial velocity 

And, finally, the Poincare variables, denoted by x, are 


(II 1-2) 


[:] 


L 

*1 

to 


(II 1-3) 


with x = ^ + g + h L = L 

= [2(L - G) ] ^ cos (g + h) = [2(L - G)]^ sin (g + h) (m-4) 

r\ 2 = [2(G - H) ] ^ cos h | = [2 (G - H) ] ^ sin h 

Note that equations (III-4) give the transformation from Delaunay to 
Poincare / , and that no singularities are introduced in this transformation. 
The inverse transformation is given by 
-1 ^1 

H = \ - tan — L = L 


-i 

g = tan — 


h - tan 


-i h 


, r 2 2 

i -1 i 2 6 1 + \ 

t&n _ G = L - 

,2 2 

^2 + ^2 

H = G - -=— — - 


(III-5) 
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In the transformation (III-5), the equations for the momenta L, G, H 
exhibit no singularities, whereas in the coordinates i,g,h singularities 
* will arise for zero eccentricity (t^ = 0) and for zero inclination 

(t) 2 = 0). This feature of the transformation will be important in later 
analysis . 

If one denotes a general canonical set of variables by w, then 
the equations of motion take the form (see Eq. (II-6)): 

• d 

w = $> — 3C(w,€) (III-6) 

° dw 


"a “I a 

.P J P 


with w = r 0C “1 0i the generalized coordinates 

the associated momenta 
where for the artificial Earth satellite problem the Hamiltonian is 
written as 


3C(w,e) £ + € 3C (1) (w) + e 2 3C (2> (w) (III-7) 

21* (w) 

The oblateness coefficient J has been taken as the small parameter e 

for convenience. Since all of the higher harmonics in the expansion for 

2 

the Earth's field are of at least one can re P resen ^ their 

2 

contribution as an e term (Eq. (II-4)). 

Now, apply a stationary canonical transformation to define a new 
set of variables w*: 


( 1 ) 


( 2 ) 


a = a* - e (3\a) - e D- t O f ,a) 

p p 


(l) 

' a 


2 (2> 

(p’ ,a) + e' E>~ (p',a) 


(III-8) 


, (i) 


which has been truncated with the second order terms. The D are the 

"generating functions" of the canonical transformation. The new 
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Hamiltonian is then 


X 


2 / (1) (1) 

'(w\e) =K(w,e) = - ^ + e K<1) + ( V D a* “ V D 

2i r l \ 


“) 


/ (1) (1) (1) (1>\ 2 / (2) (2) \ 

^ 2) + (%• fl a' ' *a' V ) + ^3 (V D 5’ " L «’ P* ) 

(1) \ 2 2 / (1) O) to) 

,«g. ) + >(- V D jS* +L a' B p'«'Vj 

’ ) /^’P’ 

W' p* 


q 2 / (1) 

_ 3tL. ( L . D~. - L 

4 
2L 


4 \ p' a 


, 2/(1) (1) 
+ 1 H_ f D - D 
2 3 X a* * p 



\U .w„> 

a’ 1 f 


(III-9) 


All functions in Eq. (III-9) are to be evaluated at w* . Choose 
D (1) (a,p’) and £> (2) fo,p’) so that 3C'(w') contains no short period 
terms except in f(w',e). This requirement is defined by 


[ 3C ' (w' , e) - e 3 f (w ' , e) ] = 0 


(III-10) 


with £’ the Delaunay variable conjugate to L' = L(w') . The Poisson 
bracket 

[A,S] = A^, - A a ,B~, = A w ,*B~, (HI-11) 

is easily shown to be invariant under a canonical transformation. In 
particular 

. _ n ^ 

(III- 12) 




then if one writes 


X 


(i) (w') = R (i) (w>) + ft (i) (w') 


with 3C (1) (w’) = av 3C (1) Cw’) 
£ * 
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one chooses 


23 a) (w’) 


such that 


J d» (1) 

?3 —Sir- 

Li 


- JC 


( 1 ) 


(HI-13) 


This defines D^(w') uniquely up to an additive function of the 
Delaunay variables other than i x . It is then convenient to choose 
to be identical with Brouwer's (L' , G' , H, £,g , — ) expressed 

as S^Ca^'). Note that the function is non-singular for zero 

eccentricity and/or inclination and is a function (as Brouwer writes it) 
of both L,G explicitly and implicitly through e and f, the true 
anomaly. When computing the required partial derivatives for £ and g 
short period variations, the singularity for zero eccentricity, for 
example, arises in the following way 


( 1 ) 


£>s 

"5U 


as 


(1) 








(III-14) 


As shown in Ref. 10, no ■— r terms arise in the case of the Hill 

e 

variables; however, zero inclination singularities still exist. That 
no singularities occur for the Poincare^ variables can readily be 
demonstrated. The argument is given for the variable \ ; similar 
arguments apply to the other variables. The function is given 

explicitly as (e* , f * , g* ,0' ,H' ) so depends on L' also 

through e* and f ' . According to the von Zeipel procedure the first 
order short-period variations of \ are given by 
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(I II- 15) 


which then can be written as (dropping the primes for convenience) 


^1 “5 l 


de ds (1) ac as (1) sh as (1) 

+ §I _ SS — SI a g SI 5h 



_ n(u - 1)(ti + 1) _ - eT \ (I II-17b) 

eLCrj + 1) 

The derivatives with respect to e and f explicitly introduce no 
singularities and neither do the last two terms of Eq. (III-17a) . Thus 
is well-behaved as e(and i) 0 . 

Next, one can show that the second order long -period Hamiltonian is 

independent of the particular canonical variables used and, furthermore, 

( 2 ) 

that the second order generating function 2) is non-singular for 

( 2 ) 

zero eccentricity and/or inclination. Recall that the function £> 
is chosen so as to cancel all second order short-period terms of the 
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Hamiltonian. In order to obtain the desired results, note that 


a) (i) (i) a) 

K - X 

p' a' a 


“ s':’ .«>] 


is an invariant for canonical variables. Furthermore, 

<2) (2) . d (2) 

V D a' " L a’ D p’ = "SI 5 


and 




(III-18) 


(III-19) 


(111-20) 


One can rewrite (from Eq. (II I- 9)) 

(1) (1) 


(1) (1) 


" V S S'a’ s p' + L a' S p'a' S p’ 

as follows (where a are understood to be the primed variables) 

(1) U) (i) . (l) (i) (l) (i) 

"V3a S p + Vpa s p ="2 Yea s p + J Vpa s p 


■ (■• i--i) (“■"').! 


(l) (l) x (i) (i) 

L p S ap S a "2 Vpp s a 


5 

/ (i) 

(1)\ 

(1) , 

£ 

/ (1) 

u>\ 

53 

is. - 
y P a 

L s~ ] 

a p J 

' + 2 

5p 

L S~ 

\ P a 

■ L a S p / 


+ 2 ( ^ L pa ~ L Qa 


( 1 ) 


\ (l) / a) ( 1 ) \ C 

I S~ - I s- L~ - S-. Lv S- 

7 P 2 \ a PP P Op/ a 


(III-21) 


( 1 ) 


( 1 ) 


/ (D (1) \ (1) / \jv vj./ t‘ v. 

2 \ S cc % • s p ^7 % " i ( S a % 


(l) 
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then also 



■r)/* 


L~ 

pa 


(l) 


I~ 

ap 


oa 


(i) 

- s~ 



Thus it can be seen that Eq. (III-2! 
Eq. (III-2I). The second order part 


(III-22) 



[) cancels the last two terras of 
of the Hamiltonian (III-9) 


consequently is given by 


2 

e 






2 ( 2 ) 

The only term of Eq. (III-23) apart from *> that depends on the 

L' 3 b J ( 1 ) ( 1 )\ 

particular canonical variables used is - £ I , which is 

( 2 ) 

necessarily short period, and £) , of course, is chosen so that the 

2 ( 2 ) 


term D , cancels all short-period terms. 

L * 3 £ 

From this invariance property of the terms in Eq. (III-23) one 
deduces that the difference between the second order generating 
functions £)^ 2 ^ of two different sets of canonical variables, such as 
any arbitrary set w and the Delaunay set y for example, -will be given 

(2) (2) t* 3 / (1) (1)\ Ti 3 / (1) (1)\ 

£) — S v ^ = _ — f s , S~, J Is S~ J + arbitrary long period term 

2n 2 \“ P / 2ji 2 \ Q ' (III-24) 
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Equation (III-24) gives a convenient algorithm for computing the 

( 2 ) 

generating function D for any set of canonical variables. In 

order to assure that all of the functions arising in the error bound 

(2 ) 

determination remain bounded, one must establish that D contains 
no singularities. This may be done utilizing the known results of 
Izsak (Ref . 9) , Brouwer (Ref . 3) and Kozai (Ref . 4) . 

In Kozai’ s paper, the expression given for the component of 

( 2 ) 

the function S shows the factor 1/e for the trigonometric 
arguments sin f, sin (f + 2g) , sin (3f + 2g) , sin (3f + 4g) and 
sin (5f + 4g) . The appearance of this factor is unnecessary and a 
suitable rearrangement of terms eliminates it. Such rearrangement 
will be shown explicitly here for the coefficient of sin f; the other 
terms can be treated similarly. The coefficient of sin f as given by 
Kozai is (omitting a nonsingular multiplying factor): 

^ £ 9(11 - 30 e 2 + 27 e 4 ) - 8 n 2 (i 7 - 38 e 2 + 11 e 4 > - 4 n 3 o - 3 e 2 ) 2 x 

(3 + T) 2 ) + t) 4 <53 - 130 e 2 - 11 e 4 )J (111-25) 

, „ H 

where 9 - — = cos i 

G 


Equation (III-25) can be rewritten as 

~ [99 - 270 Q 2 + 243 e 4 - 136 + 304 Q 2 - 88 0 4 - 4t) 3 (1 - 3 0 2 ) 2 (3 + T) 2 ) 
+ 53 - 130 Q 2 -11 G 4 - 2e 2 (121 - 282 Q 2 + 33 0 4 ) + e 4 (53 - 13 Q 2 - 11 


2 4 

or dropping the e and e terms and combining: 

2 2 

4(1 - 3 e ) r . „ 3 5-, 

[4 - 3t] - T] ] 
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Now, this can be rewritten as follows: 


4 - 3r) 
e 


2 


5 

X- 


^ J ^ h 4 + Tj 3 + 4r) 2 + 4 t) + 4] 


= (1 - 1 l )(1 + a) [r) 4 + ,3 + 4T, 2 + 41) + 4] 
e (1 + ri) 


= j— - [tj 4 + tj 3 + 4tj 2 + 4t] + 4] (III-26) 

which remains bounded as e -» 0 . In a similar manner the other 

expressions given by Kozai (for higher order harmonics J^) may be 

( 2 ) 7 

rearranged and thus it can be shown that S contains no 1/e factors. 

Of the terms on the right-hand side of Eq. (III-24) , the first is 
known to be bounded (III— 17 a,b); the second can be shown to be bounded 
by the above technique of rearranging. 

The (new) canonical variables w' satisfy the differential 
equations 

w' = 0 (w f ,e) (III-27) 

o w 

where one can write the Hamiltonian in the form 


2 

3C'(w',e) = ^ + e 3C (1) ( W ’) +e 2 K (2) (w') +e 3 s(w’,e) (III-28) 

2L ,2 (w') 

an analytic function of the variables w' and the small parameter e . 
Define next a transformation (canonical) to the "secular" variables 
w H by the truncated expressions 


a" = a* - e s-.,^",^) 

p 

P" = p* + e s~ f (p ,, ,a l ) 


(III-29) 
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where* S* is chosen so as to cancel the long-period part of the 

Hamiltonian (except near critical inclination) . Then S* and st 

p a* 

, give the first order long-period variations of a and £ i.e.: 


-L 

2« 



0 




(HI-30) 


The governing differential equations for w" become then 

»” = I X~» (w",6) + e 3 

and the solution can be written in the form 


* 



(i) 



„ <2) 



- eS^„ 

(p n 

,a* ) 

- G V 

<p* 

,a) 

- 

(p* 

,a) 




(1) 



„ <2) 



+ eS* t 

a 1 

<p" 

,a*) 

+ eS~ 
a 

(p* 

,a) 

+ € 2 ®~ 
a 

<p* 

,Q!) 


where of = a” - € S*„ (p",a' ) 

p’ = p” + €sJ, (p",a>) 


(IIX-31) 


(III- 32) 


In the definitions of what constituted long-period and/or short-period 
variations, the Delaunay variables were used explicitly (see Eqs. (III-30) 
and (III- 10)). If the von Zeipel technique is carried out for the 
Delaunay variables, then it is found that P"( = p") are constants, the 
"secular" Hamiltonian is a function of P" only and the coordinates 
Q"(= a") have constant rates. If one is dealing in any other canonical 
variables, for example the Poincare'' set x, then x" are defined to be 
the same functions of y" as x are of y. 


The function S* can be chosen to be identical to Brouwer’s long- 
period generating function considered as a function of w (see Ref. 10). 
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In this section it has been established that the generating functions of 
transformations (III-8) and (III-29) and their partial derivatives are 
bounded. Note that Eqs. (III-8) constitute transcendental equations for 
w* which may be written in the form 


a’ = a + e ^(p 1 ,a) 
p' = p + € £ 2 (p' ,a) 

From the general bounds on w one has 

|p| < B 


(III-33) 


(III-34) 


The functions ^(p' ,a) depend on trigonometric functions of a . 

Suppose now that p is bounded within some region R , and specifically, 
that p is bounded away from the boundary of R by at least eA/l-eK 
Where A and K are defined by 


U 2 (p,a) | < A 

and |£ 2 (P 1 ,a) - < K |P 1 - P j I for all , pj 


(III-35) 
in R 


Assume further that eK < 1 ; this in effect imposes a restriction on 
how close the energy and angular momentum may be to zero. Now assume the 
following iterative algorithm for computing the primed variables p* : 


with 


then 


V'n + 1 = P + 6 £ 2 ( ^> a) 

a ! A 

- P 


(I 11-36) 


lei - el < ea 

[ P 2 - ejJ < eK|pi - p| , since is also in R 

1^3 - e^l < <*le^ - Pil • since is also in R 
|p; : p;.J * eK 'p;-i - p^ 2 i 


(in-37) 
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and Ip; - pI < Ip; - p;l + Ip; - p| < |p; - p|a + eK ) 

l?3 - pI < Ip; - p;| + |p; - pI < |p; - p|d + eK + < £ k) 2 ) 


n-1 


Ip; - pI < Ip; - pI ^ <€K) J 

j=0 

n-1 

Ip; - el < eA ^ (gk> j 


(I 11-38) 


j=0 


Taking the limit 


iim |p; - p| < eA litn 
n oo n -> oo 


I 


J=0 


(eK) J 


lim |p; - p| < eA/l-eK 

n — > oo 


(III-39) 


A similar argument can be applied to the long-period transformation (HI- 
29) to deduce that p" will remain (sufficiently) close to p» . As a 
consequence of the above, one has 

|a f - oc\ = e|£ 1 (£ , ,a)| < e^ 1 (III-40) 

Also, note that the iterative procedure converges, i.e. 


i*; + i - p;i ^ £A(eK)I 


so that as n oo , | p | ^ 0 if eK < 1 . 


(III-41) 
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IV. ERROR BOUNDS FOR THE AXISYMMETRIC PROBLEM 

All of the information necessary to derive formal error bounds has 
been given in Sections II and III. One can proceed then in a straight- 
forward manner to derive the bounds utilizing known theorems from the 
theory of differential equations. However, it turns out that because of 
the nature of the differential equations, the bounds obtainable in this 
manner prove to be unsatisfactory for time intervals of the order of 
1/e . This fact is a natural consequence of the existence of a rapidly 
rotating phase in the governing system of differential equations; 
however, since only one such phase appears in the case of satellite motion, 
one can circumvent the difficulty by appealing to a known integral of the 
motion. In this section the conventional method of error analysis will 
be presented first and then the extension to the large time intervals 
will be given for the problem with an axisymmetric potential. 

A. BOUNDS FOR SMALL TIME INTERVALS 

In order to simplify the following presentation, some new notation 
will be introduced at this point. If A and B denote n-dimensional 
vectors then |A-b| will denote the matrix of absolute values of the 
component differences of A and B i.e.: 
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The governing differential equations for the '’secular*' variables were 
given as (Eq. (III-31)) 

W " = ^o [3C ,, (w" f €) + (w" f g) ] (IV-2) 

dw* 

from which the approximate state vector w^ is defined by 


W A = ® 0 ~ r( V £ > 

A o A 


w^(0) = w" (0) 


(IV-3) 


For convenience, Eqs . (IV-2) and (IV-3) can be rewritten in the form 
w” =A(w m ,€) + £ 3 S'(w", e) 

"A =A<*v e) 

Since the functions — X" and -A- X" satisfy a Lipschitz 
dw" 


(IV-4) 


dw" 


condition on the domain of definition of w(t), it follows that 


I A (w", e) - A(w”,e)| < k|w" - w^| = km" (IV-5) 

where k is an n X n matrix if m", the matrix of absolute values of 
the component differences of w" - w^, is n X 1 . The particular form 
of Eq. (IV-5) was chosen since a vector function, say A (w”), satisfies 
a Lipschitz condition on w" if and only if each of its components 
A i (w",t) does. Since the constants may be different, the use of the 
matrix of Lipschitz constants k can afford a more precise bound than 
that usually provided by the norm ||w” - w”I| . 

As a consequence of (IV-5) one can immediately write 

m” < km” + e 3 ¥(w",e) (IV-6) 
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where from a priori bounds on w(t) 

|»(w",€)| < W(c) (IV-7) 

Hence 

m" - km" < € 3 W (IV-8) 

dt 

which is readily integrated to give 


3 “1 

n" < ra M (0) exp kt + e Wk [exp k t - I ] 


0 < t < T , t o = 0 


However, since it was assumed that w^(0) = w " (°>» the initial error 
m M (0) = 0 and 

m " = | w" - w^J < e 3 Wk -1 [exp kt - I] (IV- 10) 


At this point, several difficulties of (IV-10) can be pointed out. 

The bounds (IV-10) prove to be unsatisfactory for Delaunay variables for 

small eccentricity and/or inclination since k contains the factors — 

and — l — If w i s taken to be the Hill set z, the zero eccentric- 

sin i 

ity difficulty is removed. Although the zero inclination singularity 
remains, for many purposes the Hill variables are a convenient set to use 
due to the relatively simple expressions for the periodic variations of 
the in-plane coordinates (see Vagners, Ref. 10). Taking w to be the 
Poincare / set x, satisfactory behavior is assured for both zero eccentric- 
ity and inclination. A much more serious difficulty occurs if one wishes 
to examine the bounds for time intervals of the order l/e . Expansion 
of (IV-10) yields, for "small*' time intervals 

°° k^ 1 t J 

cV 1 [exp kt - X] = e 3 Wt + e 3 W £ = — Jt < IV - U) 

j=2 
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However, for time intervals of order 1/e , the bound (IV-10) becomes 
very large i.e., behaves like exp 1/e . 

Assuming now that m" is at most 0(e) (i.e., bound (IV-10) is 
satisfactory) then one can complete the analysis by including the periodic 
terms. If this is done, the total approximate solution of interest here 
is written as 

W c = W A + (IV-12) 

with r(w^) giving the first order periodic parts of w as defined by 
eqs. (III-8) and (III-29) with the generating functions considered as 
functions of the double primed variables. Equations (III-32) can be 
written in the form 

w = w" + er(w") + e 2 5(w",e) (IV-13) 

then 

|w - wj = |w" + er( w") + e 2 £(w”,€) - - er(w^) | < 

(IV- 14) 

Iw" - + e|r(w”) - r(w^)| + e 2 U(w".€)| 

Since e|r(w") - V(w^)| gives the error in the first order periodic 
terms of the solution and m" is 0(e), then the term contributes error 
of second order. Thus the effect of the last two terms of (IV-14) can be 

o 

combined into one second order term e Z to account for all periodic 
errors of the solution. The error bound for ,, small ,, time intervals, 
assuming exact initial conditions, assumes the form 


l w “ * c l < € 3 Wk 1 [exp kt - I] + e 2 Z 

or, effectively, 

e 3 Wt + e 2 Z 


(IV-15) 


(IV-16) 
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The difficulty of the above bounds for t ~ l/e is a direct 

consequence of the existence of a rapidly rotating phase in the dynamical 

system. In the treatment of systems with rapidly rotating phases by the 

method of averaging, the governing equations for these phases are 

considered separately. The general result obtained then is that the error 

is 0(e) for t ~ l/e rather than 0(e 2 ) as one would expect from the 

2 3. 

truncations performed i.e., truncation of 0(e) periodic and 0(e) 
secular terms. In the following, such separation will be effected and, 
by appealing to known integrals, the bounds will be derived for all 
2 

variables to 0(e) for t ~ l/e . 

B. EXTENDED TIME ERROR BOUNDS 

The following analysis will be carried out for the Poincare 7 variables 

explicitly utilizing known results for the Delaunay variables and their 

rates. The secular Hamiltonian was defined from the von Zeipel procedure 

as being a function of the Delaunay momenta P" only (to second order) , 

hence in the Poincare 7 variables one writes 

3c" (x " € ) = - jL + e % (1) (x") + e 2 <* <2) (x"> - e 3 q>(x",e) (IV-17) 
2L" 2 

where K U) = ^ (a - 1 \ , H and G functions of x" 

4L" G" \ G / 

(Eq. (III-5)) and (H is F* of Brouwer considered as a function 

of x" . (Explicit expressions for <R (2) in terms of Hill variables 
for any J n may be found in Ref. 10, which could then be transformed to 
Poincare 7 variables if necessary.) Equation (IV-17) can be rewritten more 
conveniently as 

X” ( x ", e) = - -*£r + e ^(x",e) - e 3 q)(x",€) (IV-18) 

OT *> 2 
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where 


£Kx",€) } f + if, n f + q 2 ,^j 


The equations for the Poincare variable rates become 


= -^3 + e ^[n fc7(x",e) - e 2 cp(x",€)] 

L 


L” = e 3 (x” , e) 


3C (1) and « <2) 


contain X and 


x" = e<t fc7(x",e) - € 2 cp(x",e)] 

R o b ~„ 


with 


, $> a 4 X 4 matrix. 


The approximate variables x” are defined by Eqs. (IV-20) and 
with cp(x M , e) set equal to zero and x^(0) - x M (0). Consider 
differential equations for q” and : 

n" = c JL- <d( x " e ) - e 3 ^L 

n i £ di" ,€) 6 dq 


*'i = ■ € s^ (x "- e) + €3 f§ 

Recall that ^(x' f ,e) is given by Eq. (IV-19) so that with 

A „ d& 


N 1 = 2 


n / m 2 f t2 ^ 

0(^1 + ^ ) 


(IV-19) 


(IV-20) 


do not 


(IV-21) 


(IV-21) 
first the 


(IV-22) 
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a bounded quantity, one obtains 


< = e Vi ’ " 3 ^ 9 
1 


3 fc 


(IV-23) 


Hence 


1 d_ 

2 dt 


€N i\ + £ ^ 


t m2 m 2 3 tt dcp m dcP 

+ \ I = e 5^ ' \ 3|* 


(IV— 24) 


The approximate solution x” of interest is given by cp = 0 . thus from 
the boundedness of cp (and its partial derivatives) 

I {< - ^lA 2 ) + - *1A 2 ) I = I A + < ) I < (IV " 25 > 


Here, as well as in the following discussion, the extended time interval 
will be taken as t ~ 1/e so that 



(IV-26) 


The reader may prefer to think of the time interval as defined by 
nt ~ 1/e where n is taken to be the (suitable) mean motion. For 
mathematical convenience, the definition t ~ l/e will be used. 

Now, rewrite Eqs. (IV-23) as a single complex equation (j = \/-l) : 

(|” + Jt]”> + 

and the approximate equations as 


dcp 

K ' 


dcp 

K 


(IV-27) 


e" + i 

S 1A J 


\a = j 


n ia<*ia + 


^1A> 


(IV-28) 


N 1 = N 1 < + < • ^ + * L "' £) 


148 


ERROR BOUNDS ON POSITION AND VELOCITY 
Difference Eqs. (IV-27) and (IV-28) to get 

d'; * id'; = + jn”> ♦ + j*;> + * 3 [f^ - j fjj 

(IV-29) 

From Eq. (IV-26) : 

**1 = K - n i A i < vffi 2 + v’ 2 ) < 

and (IV-29) thus becomes, with the aid of an integrating factor, 

I £ coat; + < 

r (IV-30) 

eV^iA* ^ - J sfj + JVi ■ Vi j 1 

Then, since the right hand side of (IV-30) is bounded, it follows that 

| (£6^ + J£Jlpe _JeN lA t |< e 3 M 3 t = e 2 M 3 t ~ 1/e (IV-31) 

but | e"** eN lA t | = 1 so 

+ j An" 1 < e 2 M 3 4 (IV-32) 

From similar arguments, it follows that for and 

IAI 2 + - * K 4 • t ~ 1/e (IV— 33) 

Also, from the differential Eq.,(IV-20) 

|l" - L”| = |e 3 J ^ dt| < € 3 M 5 t = € 2 M 5 (IV-34) 

The remaining coordinate \ v causes some difficulty, since with 
| L n - L” | known to O(e^), a straightforward analysis of the 
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equation gives lx” - \^\ only to 0(e) for t - 1/e . In order to 

obtain bounds for x" consistent with those of the other coordinates, 

one must appeal to the knowledge of an exact integral for the axi- 

symmetric problem. In effect, one can re-define the mean motion as 

introduced by Brouwer (Ref. 3), who wrote 

2 


A 

n = 


(IV— 35) 


and 


hence, with J i and JL functions of L' ,G",H only 

1 4 


i" = n [l + <=JZ. + e 2 i ] + 0(e 3 ) 
o l & 


(IV-36) 


Recall that x” was defined by Eq. (IV-20) , which in a more explicit 
form is given by 


X = ' 


where 


1 - e 


3 . .n*l 3 » 2r I TJh"\ 2 , H-l 

2 k i L - e L 5 v?) - 1 ' 2 


g 3 dg 

(IV- 3 7) 


+ e “ 6 2 | “ e ” Sl" - 


2_2 

k .H 5 a 

1 2G- 3 


(■W-O 


with G” = G ,f (x n ) and H” = H M (x") 


L' 3 


( 2 ) 


5^ = — — -n — , a bounded quantity 


Define now a new constant a 


by ^with (ft* = 2L 2 "’ ^ ^ 

2 

-tL = - 3C = - 3C" = -H_ [1 _ tk l"" 1 - e 2 ft*+ eV ] (IV-38) 

2a „ r ..2 1 1 


and a new "mean motion by 
3/2 


J . (!) - 4,[i - - . 2 «* . «\] M UV- 


39) 
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or expanding 
2 

^ u 

n = 


[ , 3 . T „-l 2 3 (O* 2 3 .2-H-2 3 x 

1 - G 2 k l L “ 6 2 + £ 8 k l L + € V x ' €> 


(IV-40) 


Thus 

X" = ft - 6 ^4 


K** [„(! 

;) 2 - 1 - 2 jd 

. e L + | «* - | kV' 2 ] 

4G” 4 [ \ C 

■ J G 

[_ 2 2 8 1 J 




(IV-41) 


Again, the approximate is defined by (IV-41) with si*"" and 

cp* = 0 . From the exact known integral, H" = H in Eq. (IV-41) and from 
the definition of G": 


e ti2 m2 

It 

G" = L" - — — 


then from Eqs. (IV-32) and (IV-34) 


so that finally 


|g" - G^| < e 2 M 6 t ~ 1/e 
lx" - \”\ < t ~ l/e 


(IV- 42) 


(IV-43) 


(IV-44) 


If one is interested in orbits with non-zero eccentricity and/or 
inclination (i.e., e » e, sin i » e) then the above analysis can be 
carried out analogously for the Delaunay and/or Hill variables* In 
particular, for the Delaunay variables, the rapid phase is l" and an 
equation similar to (IV-41) (but somewhat simpler in form) results for V' 
Due to choice of g and h as the other two coordinates, the € term 
of (IV-41) is found to disappear. (Of course, the functions ant * 

q> are different than for the Poincare / variables) . 
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C. THE INITIALIZATION PROBLEM 

At this point the relevance of the above results to the so-called 
initialization problem may be noted. The two primary uses of an analytic 
(artificial satellite) orbit theory are orbit determination by fitting 


to observational data and orbit prediction from some initial state vector. 

In the case of orbit determination, the mean (double primed) variables 

are obtained to high accuracy by fitting to observational data. This 

accuracy depends on the number and quality of the data points. In this 

application, the question of initial value errors does not arise. 

The initialization problem may be defined as follows: given some 

initial radius and velocity vectors, compute a satellite ephemeris for 

some extended time interval via an analytic theory. The initial radius 

and velocity, and hence the instantaneous elements, are assumed to be 

known exactly. Analytic theories are usually formulated so that certain 

constants of the solution are mean elements, for example L' ,G M and H 

in the Brouwer theory, instead of initial values. Thus from the known 

set of instantaneous elements, the mean elements must be formed by 

subtracting out the periodic variations. Since one is considering a first 

2 

order theory, the mean elements thus defined will be in error by 0(e ). 

It can be noted here that a numerical iteration procedure has been applied 
to the determining equations (Cain Ref. 15, Arsenault, et al Ref. 11) 
which are written as 


Q = Q" 


P = P" 


(1) 

* 

1 

p, (P' ,Q) 

+ S~„ (P ,Q' )J 

(l) 

<*) I 

% (P '- Q) 

+ s~, (P",Q')J 


(IV-45) 


152 


ERROR BOUNDS ON POSITION AND VELOCITY 


Such a procedure can, of course, only remove that second order error 

that arises from considering the S functions to be functions of the 

instantaneous elements (P,Q), but still cannot account for the truncated 

second order terms. Thus from an accuracy point of view, such iteration 

procedures are of dubious value, since as shown by Eq. (IV-37) (with e 3 

terms truncated) the error in , or equivalently , will still 

2 

grow as e t from the zero order term. The other variables of either 

x" or y" do not present any problem since their rates are either zero 

o 

or multiples of e, so that an initial value error of 0(e) will grow 
3 

as e t giving results consistent with the expected accuracy of the 
trunca ted theory . 

With the algorithm suggested by the analysis of subsection IV-B, 
the initialization difficulty can be resolved. As noted, for all variables 
except the rapidly rotating phase, the use of mean elements defined by 
instantaneous value minus the periodic terms (considered as functions of 
the instantaneous elements) will lead to no difficulty. The necessary 
initialization procedure for is given by Eqs. (IV-38) , (IV-39) 

and (IV-41). The numerical value of n is known exactly from instanta- 
neous 3fC and the remaining terms of (IV-41) have at least an e multiplier. 
For the Delaunay variables, one rewrites Eq. (IV-36) with n ,r the mean 
motion defined by 

2 

n" = i" = [l - e | + *\] + 0 (e 3 ) (IV-46) 

so that, with use of energy^ 

t , 3 

The explicit expression for is identical to — — F** . [See 

, 2L’3 _ _ M- 

Vagners (Ref. 10) where it is given as — ] ]. 
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n" = £ = n + e 2 — ["■— “ k ? L * 2 1 + £ 3 terms (IV-47) 

L' 3 L 2 2 81 J t = 0 

2 

in which in the e terms one replaces the double primed variables with 
the instantaneous elements. The new mean motion n is again given by 
(IV-38) and (IV-39) . All terms appearing in the brackets of (IV-47) are 
functions already known from the general theory. 
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V. THE ASYMMETRIC POTENTIAL FIEU 
If one includes the longitude dependent terms (tesseral harmonics) 

* 

in the gravitational potential, some modifications to the analysis of 
Section IV are necessary. The additional terras in the Hamiltonian are 

00 n / R \ n 

** = r Z Z J n,m \~7 ) P„ (sin P) cos ra( *- " < V -D 

n =2 m =0 

where J \ are constants with J - 0(J 2 ) 
n » m n,m n,m 

A. = h + tan 1 (cos i tan u) - u t 

© 

W 0 the angular velocity of the Earth 

(Time is measured from an instant when the right ascension of Greenwich 
is zero.) 

In this discussion, X T will be considered first as a function of the 
Delaunay variables JC T (L,G,H,f,g,h - w # t) . To remove the explicit time 
dependence, define a new canonical variable as h* = h - w^t conjugate 
to H with the associated Hamiltonian given by 

K = 3C - W 0 H ( V- 2 ) 

where 3C now is the original Hamiltonian including both zonal and tesseral 
harmonic effects. Since time is not present explicitly in K, it is a 
constant of the motion. 

Following the von Zeipel procedure, "remove" all of the periodic parts 
of the extended Hamiltonian K via a suitable generating function, 
defined here up to second order (since is second order in e) so 

that the new variables become 
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L = L * +6 


(!) 2 as (2) dS T 

^J- + e ^T- + ^e 


dS 


as (1) 2 *S (2) , 

& = * + e TP- + € “§17“ + bv 


( V - 3 ) 


H = H* + 


bh* 


Equations similar to those above hold for the other variables. Note that 
H now contains fluctuations but that these are of second order. As 
before, and S^ 2 ^ are chosen to cancel all zonal short-period 

terms up to, and including, second order. Thus one is left with (omitting 
the € 3 function for the time being) 

2 


K ' = 


2 L ' 


+ € JCj + € (ft 


2 dS_ dS 

2 “ <2) -“®H- + ^3 < V - 4 > 


where 


3^ = 3C + 3£ t = ^ ^ S\k Ja. 6 . 1 ) 008 ^ + k 2 g + mh * + phaSe) 

( V - 5 ) 


k 2 k l 


in which 3(L includes all terms with k 1 t 0, the short-period part of 
3<^ , and 3C T gives the daily fluctuations from h* , ^ = 0 . Choose 


so that 
T 


2 Ss as 

■ + w — r = - 3C„ 


L* 3 5T 


dh 


( V - 6 ) 


Strictly speaking, (V— 6) should be written as 

2 as as as* 

— - H’ i ■ - + u> + (oj - fi ) - 

L ,3 df ® dh * ® dh * 

- - 3$j, (L ,, ,G ,, > H , ^g',h* , ) 

however, (V-6) is accurate to 0(e 2 ) since fl is 0(e). Solving by 


( V - 7 ) 
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inspection, with n* = , 


m A, 


V' V k l k 2 m *' 

T = Z Z 2 k n'-mo) sin( V' + k 2 e ' + "* + P hase > (V-8) 

k 2 k l m ® 

In Eq. (V-8) (as well as (V-6)) one can use the primed variables or 

3 

instantaneous variables with 0<€ ) error. Clearly, there exist orbits 
for which n' is commensurable with U) and hence a particular 
< k 1 n * ™ mto 0 ) 8°es to zero. These are the so-called tesseral resonance 
cases and will not be considered here. Introduction of S T in the above 
manner leaves 


K" = - 


2L 


,2 


, f « . c 2^ (2) lt 3 ^ 
+ e + € 01 — - e cp 


(V-9) 

3 


where the e function is now included, and is different from the e 
function of of the axisymmetric problem. 

The instantaneous elements may be written as a sum of the "secular* 1 
(double-primed) part and the periodic parts. In particular, with H" a 
constant, 


H = ** + “s.P + V 


From (V-3) and (V-6) it follows that 


(V-10) 


i H" = U) H + 


ds T 

§T 


thus, from (V-8) 

as,. 


(V-ll) 


r-^ ^ ^ k i nt t 

dT = Z Z z A k k m cos< V’ + k 2 8’ + «**’ + phase) (V- 12) 

f m 1 © 12 

2 k l 


Combining (V-12) with (V-5) one obtains 
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X XXf-^s* 

k 2 k l m \ ' 


cos(k,i' + k„g' + mh + phase) 
ml 2 


(V-13) 


_ mw + t 

- ) V V A cos(kJ' + k g* + mh + 

Z Z Z V k ! k 2 m 1 2 


phase) 


k 2 k l m 


The constancy of K and 3C n to 2nd order implies (because of (V-ll)) 
that 

^ ®rp P 

- + n' tj-j = tesseral fluctuation of 3C = / rpj dt (V-14) 


where J ^ dt denotes a specific second order approximation to the 
indefinite integral, considering only V and h* as time varying. 

This may be checked by forming 

5 t = T 7 = ^ I sln( V + k 2 g ’ + rnh *’ + phase) <V “ 15) 

k 2 k l “ 

If one assumes that, to the order necessary here, 


k-ji' + mh* ' « [k^n* + mu^Jt (V-16) 

then conclusion (V-14) follows. 

For discussion of the error bounds, the formalism of Section IV can 
be retained to a large extent. Due to the absence of secular tesseral 
terms, only the cp function of the secular Hamiltonian will change and 
hence the bound on that term will be different. It follows then, that 
the error bounds on x", with the exception of \ , are derived in 
exactly the same manner as before with different values for the constants 
M i . The derivation of an error bound on x" is not as simple as before, 
since for the asymmetrical problem the two separate integrals of energy 
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and polar component of angular momentum no longer exist. 
From the extended Hamiltonian (V-14) one finds 


X = 


1 - e | ^L'" 1 




3U 2 R 


4G 


[•W- 


(V-17) 


3 dcp* 

€ Sv 


The object again is to obtain an expression for (l - e * k L' - ^) 

I/ 3 2 1 

accurate up to, and including, second order. Since K is a constant 
of the motion, 

K = 3C - 0 ) H = K*’ = 5C’ f - to H" 


so that 

30' = X - ul (H - H") 

© 

From general theory for time dependent Hamiltonians 


d3C dX 

dt ST 


and 


3C 


■I 


dX „ 

dt = constant 


(V-18) 


(V-19) 


(V-20) 


(V-21) 


The constant is related to quantity 3C" , which is also a constant to 
second order as defined by the von Zeipel procedure. In fact, if one 
chooses a particular approximate evaluation of the indefinite integral, 
as in (V-14), then 

dX 


K" = 3€ - r 

J A 


St 


dt 


(V-22) 


Note that as defined by (V-19) , 3C M is known only to second order with 
third order secular terms (from H"). After times of order l/e this 
constitutes an error in 50** of second order secular and thus finally to 
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first order in if L' is defined through 3C" . Using the relation 

(V-22) avoids this difficulty, since it turns out that the third order 

evaluation of the specific integral / &K/dt dt yields terms that are 

J A 

still third order after t ~ l/e, with tesseral resonance situations still 
ruled out. This may be verified by considering the first order variations 
of the variables in cWdt and noting that m ^ 0 in (V-15) . 

Next, define as before (with 3C" defined by (V-22)) 


aS = 


3C" 


and the mean motion by 


n = [p/a 3 ]^ 


(V-23) 


(V-24) 


so that 




£ I k i L 




2 3 , 2 T , - 2 2 3 

+ 6 8 k l L ‘ £ 


§ <ft*j + e 3 tp* (x",e) (v-25) 


and finally, the equation 


,-3 






(V-26) 


The approximate solution is again defined by (V-26) with the Cp* 

functions equal to zero. The argument for obtaining the error estimate 
follows as before. 

The algorithm for computing the correct initial value of the mean 


r ws t 


dt 


motion n now involves the evaluation of the integral 

A 

This may be done by a suitable expansion on eccentricity; one such 
evaluation is given in the Appendix. 
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VI. NUMERICAL VERIFICATION 


Equation (V-26) then provides an algorithm for computing the correct 
initial conditions (to the order of accuracy demanded by the general 
solution) in the case when the total potential of the Earth is taken into 
account. The algorithm includes the (suitable) evaluation of the 

r ^ 

indefinite integral J dt . It is of interest to obtain numerical 

verification of the general accuracy theory of Section V . The explicit 
p ^ T ^ 

expression for J dt has been derived earlier by Vagners (Ref. 16), 

and was subsequently incorporated into the Lockheed Closed Form Orbit 

Determination Program (Ref. 17). This program utilizes a complete first 

order analytic solution that is equivalent to the extended Brouwer solution. 

(The extended Brouwer solution is taken to include J 2 short-period, 

J? and general J long-period, J medium period (daily) effects and 
^ Xi n, m 

all second order secular effects not accounting for tesseral resonances 
(c.f. Giacaglia (Ref. 18), and Garfinkel (Ref. 19).) The Lockheed solution 
is due to Small (Ref. 12), and Vagners (Ref. 16). 

Since the error in the mean anomaly (or equivalently \") is directly 
related to in-track position error, the simplest test of overall accuracy 
is to compare the in-track, cross-track and radial positions as predicted 
by the analytic solution and numerical integration of coordinates. Since 
the time intervals of interest are of the order 1/e , the comparison was 
performed over a seven day interval. The test orbit was of 2000 nautical 
mile altitude and circular. Such an orbit, not including the results 
of Section V (roughly a 200 foot error in the semi-major axis), resulted 
in a 200 mile in- track error*. After "tuning" the mean motion with the 


t The comparison study was carried out to determine the effects of inclu- 
sion (or omission) of tesseral harmonic short period terms in the semi- 
major axis. The energy had already been incorporated in the formulation of 
the axisymmetric problem. 
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energy, the secular error was decreased to 900 feet, which is an order e 

decrease as demonstrated by the theory of Section V. The comparison is 

shown in Fig. 1, where it can be seen that the periodic errors and the 

2 

secular error are now of the same order of magnitude i.e., 0(e ). The 

cross-track and radial errors are periodic and have amplitudes of ± 120 
feet and ± 350 feet, respectively, for the comparison orbit. 
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VII. CONCLUDING REMARKS 

In the previous sections error bounds were derived specifically for 
the Brouwer procedure using the Poincare / variables. From the general 
theory, an algorithm was derived for the correct computation of the 
initial conditions for the Brouwer theory. It is then of interest to 
note the relevance of the results of this paper to other orbit theories, 
and also to present the computation of coordinates from the Poincare 7 
elements . 

Insofar as the first item is concerned, exactly equivalent errors are 
to be expected from any complete first order theory provided that care is 
taken in establishing the correct mean elements for that theory. A 
complete first order theory is defined as one that includes the first 
order periodic and second order secular influences of any harmonic. This 
distinction is necessary if one wishes to compare theories for prediction 
of orbits from a fit to observational data or for prediction from an 
initial state vector, i.e. the initial value problem. For example, the 
theories of Kyner (Ref. 1), Petty and Breakwell (Ref. 21), including a 
time equation carried only to first order secular terms, would give 
satisfactory results if applied to orbit prediction from a fit to data. 
However, for the initial value problem, these theories would prove 
unsatisfactory (giving e errors for time t ~ 1/e). The latter 
difficulty could be remedied if the time equation (or its equivalent) 
would be carried out to include second order secular effects and an 
energy algorithm used to calibrate the mean motion. The theory of Small 
(Ref. 12, Ref. 16) is a complete first order theory and includes the 
correct algorithm for computation of initial conditions. 
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With more or less difficulty, any theory appearing in the literature 
may be analyzed in a manner analogous to that given in this paper and 
equivalent results obtained. In each case, the energy will have to be 
used to establish the mean motion (or the constant rate of the fast 
variable) to second order, unless complete second order periodic 
expressions for the semi-major axis are available. The questions of 
error bounds become more difficult if one admits orbits at critical 
inclination and/or orbits at resonance with the tesseral harmonics. Such 
orbits are excluded from the general class investigated in this report 
and remain the topic of future investigations. 

The last point to consider is the computation of the coordinates 
from the Poincare^ elements in which most of the theory of this paper was 
developed. In terms of conventional orbit elements, 


^ = M + to + ft L = y/\±a 

TJ = [2 Ji±a (1 - \I 1-e 2 ) ]^cos(to + ft) = [2 (1 - / 1-e 2 ) ]^sin(to + ft) 

T) 2 = [2 ” cos i)]^ cos ft i 2 = [2 \/\±P (1 - cos i) ]^sin ft 

2 

where M is the mean anomaly and p = a(l - e ). The remaining elements 
were defined in Section III, Eq. (Ill— 1) . Known the time t one can 
find and L » then compute 


e cos (to 


\ f £i + 

+ fl) = L i l 1 " 4L 

S -^1 


(VI 1-2) 


e sin (to + ft) 

An iterative procedure yields A, e cos f, e sin f defined by 
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tan 1 (e sin f) >J 1-e 2 (e sin f) 

A = 2 j - + (VII-3) 

n 1 +vl-e 2 + (e cos f) 1 + (e cos f) 

n n 

(e cos f) = (e cos (w + ft)) cos (\ + A ) + (e sin(w + ft)) sin (\ + A ) 
n+1 n n 

(e sin f) = (e cos(w + ft)) sin (\ + A ) - (e sin(u> + ft)) cos (\ + A ) 
n+1 n n 


where 


•J l-e^ = 


2 2 

jl + \ 

2L 


3o that the radius is given by 


L 5/2 (l - e 2 ) 


where 


D = 


l* + 


1 - 


2 2 
^1 * ^1 
4L 


i 


[\ cos 


and the cartesian coordinates x,y,z by 

3/2 2 i poo 

x = ( 2D ~ [<2L - l 2 ) cos ( X 


L 3/2 (l - „ 2 )* 
2D 


lt 2 \ COs( ^ + A) + (2L - ^ 


- - <«• - - *4 - - *>* 


(VII-4) 

+ A) + sin(x, + A)] 

+ A) + n 2 S 2 sin (x + A) ] 

- \ - n 2 ) sin (x + A)] 

[r) 2 sin(x + A) - cos (\ + A)] 


165 



TRAJECTORY ANALYSIS AND GUIDANCE THEORY 


o 



I I 


30N3H333I0 MOVdl-NI 


166 


COMPARISON OF NUMERICAL INTEGRATION 
AND ANALYTIC THEORY ORBIT PREDICTION 

FIGURE 1 

jproduced from LMSC Document No . 579774 
icv of Lockheed Missiles and Space Company 
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APPENDIX 


Explicit Evaluation of 


/. 


oK t 

“St" dt 


integral 


For the evaluation of the initial value problem, the indefinite 

r 

j A ~ s *~ dt must be evaluated or, equivalently, the generating 
function S T must be found. It will be assumed here that the integrand 
is given by (V-15) and the integration will be carried out in conventional 
variables. 

The following expressions prove useful: 

i (D <- » £ ••• “t 2£ ... 


£=0 


where 


n-m 

2 

n-m-1 


for n-m even 


for n-m odd 


m-s s s 

cos u sin u cos i 


<- ir X 


cos 0 


where 


■<* - V.’ ■ I (!) 

s=0 

[r i cos m (h* - x*) + r 2 Sin m(h* - \*)] 

s/2 if s is even; = 1 , = 0 


+ 1 


if s is odd; r i - 0, T 2 = 1 


(A-2) 


X, = a + x with a the right ascension of 
o n, m o 

Greenwich at a base time t 


169 


TRAJECTORY ANALYSIS AND GUIDANCE THEORY 


n+1 

+ - I (T) 

r p=0 


2RC1 - e Z ) 


2.p+l n+1 


t' 

hT S (q) cos ^ p “ 2c i^ u “ (A- 3) 


q=o 


n j u cos k u = J J (c) (d) Tj *T [5 1 C ° S(j + k - 2C “ 2d)U 


c=0 d=0 


where 


i = 


+ & 2 sin(j + k - 2c - 2d)u] 


j + c + j/2 if j is even; & = 1, & 2 = 0 


(A-4) 


j + c 


+ if j is odd; =0, 


6„ s 1 


r 

At this point, the assumptions under which J dt will be 

integrated may be stated. The inclination angle i and the eccentricity 
e will be taken as constants. Since no appreciable difficulty is incurred 
thereby, the following will be adopted 


with 


U) = U) + OJ* u 
o 


h* = ft + H'u - u t 

o ® 



(A- 5) 


The last item is the central angle-time relationship. Since the integrand 
is (essentially) now a function of u, one would prefer to integrate 
with respect to u. To the first approximation 

~ 2 

du = n dt + 0(ee ) (A-6) 
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so that the contribution of the J term is given bv 

n, in J 

/ n-m-2£+s m-s 

sin u cos u [tj sin m(h* - x *) - ^ C os m(h* - x *) 


] X 


cos(p - 2q) (u - u) du 


with 1 * 1 


n, m 


v n+1 p 


H J 


n.m n+1 _ 
a 2 


tr I I I 2 jSggf o (T; 

5=0 p=0 q=0 s=0 


(, p )( :) (- ■) ti! 


2 P (1 - e 2 ) P+1 


g n-m-2f, 

X cos i sin i 


(A- 


or with h = h* - x,* , 


E n,m /-CP - 2q)(u - sin nth - r 2 cos mh]^ cos(j + k - 2c - 


with 


Then let 


5 2 sin(j + k - 2c - 2d)u]du 




- 1 ) 
j+k 


c=0 d=0 


B o = - Wq (p - 2q) 

B x = (1 - w'Xp - 2q) 
B 2 = j + k - 2c - 2d 

s 3 - .<«, - 1*) 
■4— (-'♦*) 


so that the integrand becomes 


If one prefers, the F and G functions of Kaula, Ref. 22. may be 
used instead. 


> x 

■7) 

2d)u 
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(A-8) 


J cos (B^ + B 1 u)[& 1 cos B 2 u + b 2 Sin sin(B 3 + B 4 u) 

_r 2 cos < B 3 + B 4 u >l du 

The following non-zero combinations arise in the above integral: 

I = / cos (B + B u) cos B_u sin(B Q + B u) du 

1 J o 1 2 J 4 

I = - / cos (B + B u) cos B u cos (B + B u)du 

2 J o 1 2 

I = / cos (B + B u) sin B u sin(B + B u)du 

3 J o 1 2 0 4 

I = - / cos(B + B.u) sin B u cos(B + B u)du 

4 J O 1 2 2 4 

which can all be evaluated explicitly to give 
rrn- cos[B 3 -B o+ (B 2+ B 4 - Bl )u] + B 


(A-9) 


h-i 


B r B 4- B 2 


14 2 


2 4 1 




* B cos[B 3+ B 0+ (B 1+ B 4 -B 2 )u] - * cos[B 3+ B o+ (B 2+ B 4+ B 1 )u] 

~ ' 1 4 12 J 

dhr sin[B 0 -B 3+ (E 2 -B 4+Bl )u] - B -B -B 
4 1 14 2 

ZF-Ib7 sin[B 3 + B o+ CB 1 + B 4 -B 2 )u] - - V B o +<B l + VV u] j 

B -hr -n[B 3 -B o+ (B 2 + B 4 -B 1 )u] + BBB sin^-B^CB^-B^u] 

1 4 2 2 4 1 

TFTb" sin[B o+ B 3+ (B 1+ B 4 -B 2 )u] - B +B +B ~ sin[V B o +( WV u] 


B 2- B 4- B l 


X 3 " 4 


2* 4 1 4 12 


B -B -B ‘=o S [B o -B 3+ (B 2+ B 1 -B 4 )u] + B B +B cos [ B o +B 3 +( V B 4 _B 2 )u]+ 
4 1 2 4 2 1 
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+ W B 2 COs[ VV (B 2 + W U J - 5-i i+ B 2 COs[ V B o +(B 2 +B 4 +B 1 )u 3 


In the expressions of Ref. 16, it was assumed that i, u ,n and r 
were constants and in the test case the orbit was circular. The exten- 
sions of the above development cause no difficulty other than increasing 
the number of terms. However, any improvement of accuracy for non-zero 
eccentricity orbits is difficult to assess due to the approximation of 
Eq. (A-6). In order to define the error remaining as of order e 2 € B 
one must include the e terms in (A-6) . 
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An Investigation of High Eccentricity 
Orbits About Mars* 


by J. V. Breakwell and R. D. Hensley 

Stanford University 
Stanford, California 


H67-29377 

SUMMARY 

Possible long-period fluctuations in the radius of pericenter, r , of 

P 

an orbiter of Mars due to the solar gravitational field are investigated. 

The study is restricted to a "region of interest" defined by 6000 km > r > 

P 

4000 km. Eleven different "critical" orbital inclinations are found for 
which the long-period fluctuations in eccentricity and inclination contain 
terms of vanishing frequency suggesting very large amplitudes. A closer anal- 
ysis of these resonant situations near a critical inclination is accomplished 
by transforming the Hamiltonian into that of a simple pendulum problem. Maxi- 
mum variations in the eccentricity, and hence radius of pericenter, are then 
obtained and curves of maximum change in radius of pericenter versus eccen- 
tricity plotted. 


I . INTRODUCTION 

For many of the problems in space-flight mechanics it is necessary to 
find an orbit which satisfies certain boundary or mission conditions. In this 
investigation, which considers the problem of an artificial satellite about 
the planet Mars, the size and shape of the orbit are dictated by two practical 
considerations: (l) the savings in fuel obtained by transferring from an in- 

terplanetary trajectory to an elliptical orbit and (2) the desire to pass as 
close as possible to the surface of the planet. This last consideration intro- 
duces the notion of the radius of pericenter r^ which, with the eccentricity 

★NASA Grant NsG-133-61 
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e, will determine the size and shape of any orbit. The mission requirements 
can be expressed in terms of these two parameters, e.g., the orbits must lie 
in a range of interest for which e > 0.5 and 4000 km < r^ < 6000 km (see Fig. 
l). Of major concern to the mission planner is the variation of r^, over a # 
long period of time, caused by a fluctuation in e. 

This study investigates the possible long-period fluctuations in the 
radius of pericenter of an orbiter of Mars due to the solar gravitational 
field. For those orbits which satisfy the mission requirements the "secular" 
rotations of the orbital plane and the major-axis orientation (argument of 
pericenter <n) due to Mars' oblateness dominate those due to the solar field. 
On the other hand, these oblateness secular rotations are substantially slower 
than Mars' motion around the Sun. This latter consideration permits averaging 
of the perturbations not only over orbital revolutions but also over the Mar- 
tian year to obtain equations describing "long-period" rates of change of in- 
clination I, eccentricity, longitude of the node 0, and argument of peri- 
center cjj, the last two elements being driven mainly by oblateness. The 
long-period rates of inclination and eccentricity are sinusoidal in certain 
linear combinations of oo and ft leading to fluctuations in inclination and 
eccentricity as sums of easily computable sinusoidal terms of known frequency. 
There are 11 different "critical" orbital inclinations where the frequency of 
a sinusoidal term vanishes leading to a very large amplitude. 

An analysis of these resonant situations near a "critical" inclination 

I is accomplished by transforming the Hamiltonian into that of a simple pen- 
C (l) 

dulum problem. And, as in the pendulum analysis, it is possible to obtain 

a phase-space contour describing the total libration of the system which, for 

this study, determines the maximum excursions of I and e near 1^. Once 

the maximum variation in e is known the maximum change in r (&r ) can 

p p max 

be calculated, for those orbits meeting the mission conditions, and curves 

showing (6r ) versus e plotted, 

p max 

As noted above, in this study the radius of pericenter lies in the range 

4000 km < r < 6000 km 
- P 
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Introducing the semimajor axis a and the eccentricity e, a "region of 
interest" is defined by 


4000 6000 

TT^y< a <xr^y 

This region is shown on Figure 1. 

II. THE DISTURBING FUNCTION 

In this investigation the motion of a body in the gravitational field of 
the planet Mars will be studied under the assumption that the only perturbing 
forces acting are the oblateness of Mars and the solar gravitational field. 
Ignoring, of course , the effects of the (small) mass of the orbiting body, 
then the equations of motion in a Mars-centered coordinate frame may be 
written in the form 


r = V( V R g ) 


( 1 ) 


where r is the radius vector of the satellite, 7 the gradient operator 
and V R s the potentials of Mars and the Sun respectively. The solar per- 
turbing potential R may be written as 
s 


R 


S 




( 2 ) 


where [i a is the gravitational constant of the Sun and r fl the radius vec- 
tor from Mars to the Sun. The potential of Mars can be expressed as 

R=^ + R = 

M r m 

with J 2 the second 
equatorial 




(l - 3 sin 8) + higher order terms! (3) 


harmonic coefficient 
radius of Mars 
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8 the planetocentric latitude 

R the perturbing potential due to the figure of Mars 
in 

The higher order terms indicated in (3) will henceforth be omitted and 
will be taken as the perturbing potential due to Then the disturbing 

function as referred to henceforth will be given by R m + R s * Equation (3) 
then is easily rewritten in terms of the satellite's inclination I, (mea- 
sured from the Martian equator), true anomaly v and the argument of peri- 
center cjo as 


^1 - | sin 2 I + | sin 2 I cos 2(v-kjd))j (4) 

Since the effect of the terms periodic in the satellite ' s mean anomaly M ar< 
not significant, R^ can be averaged over 2jt on M. Transform, therefore, 
the true anomaly v to the mean anomaly M by the differential equation 




(5) 


and express the radius as a function of the semimajor axis a, eccentricity 
and true anomaly, 


a(l-e 2 ) 

1 + e cos v 


( 6 ) 


so that by obtaining from (5) and (6) 


the relations of Tisserand 


( 2 ) 




dM = 


(1-e ) 


1 

2,3/2 



sin 2v 



cos 2v = 


0 


(7 


the averaged disturbing potential due to J ^ becomes 
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R 

m 


2jt 



dM 


J 2 R e / 
2a 3 (l-e 2 ) 3 / 2 1 


3 4 2 \ 

2 Sln V 


( 8 ) 


Note that R is also independent of the argument of pericenter a), 
m 

Next, consider the expression for the solar perturbation potential R . 

_ — s 

If B denotes the angle between r and r g , then the first term of Eq. (2) 

can be developed in a series of Legendre polynomials as follows: 



so that, together with the fact that 

r • r 

cos B = — (lO) 

rr ' ' 

s 

the potential R g becomes 



From Fig. 2 the unit vectors e and e , in the radial direction of the 

r s 

satellite and the Sun, may be written in component form as 


cos ft cos 9 - sin ft cos I sin 9' 
sin ft cos 0 + cos ft cos I sin 0 
sin I sin 6 


and, with C = cos and S = sin, a convention used henceforth: 


C A C T 

A L 


C T S. 
L A 


( 12 ) 


(13) 
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where Cl is the right ascension of the satellite's ascending node 

0 is the satellite central angle 

A is the right ascension of the Sun 
s 

L is the plane tocentric latitude of the Sun 


Thus the Sun's disturbing function becomes 
2 


U r 


R = 


2r 


l [V“- 


■A C 8 - S n-A C L c i s e + S L s 


4 


R is next written in terms of the satellite's true anomaly via 


(14) 


so that 


where 


= C 

0 V+CJO 

S = s 

0 V-KJO 

C 8 = C v< - Wto + « 

S e C e ■ Wife + I c v s an - I c v s au 


U r 


2r 


" ■ Kv-0 2 ■ ( c ‘V°-\ ■ V ')1 

* = k S L S I C a-A - c l s c iVaVaJ 

L S S s S S SJ 

c ■ K c k * - v-)1 


(15) 


{f a S ( ™ ) + 3Ss 2(v W ) + § e - 1 } (16) 
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Holding the Sun's variables L g and A g constant and integrating R with 
respect to time, the short -period terms depending on 0 can be removed. The 
average then is 


(17) 

which may be readily integrated by introducing the eccentric anomaly E 
through 


- i r 2fl 

R ° = ^Jo *' 


dM 


dM = (l - e C ) dE 
£ 


and the relationships 


C 

v 


< C E ~ e > 

1 - e C 


(2 - e ) C - 2e 


C + 2e 2 - 
E 


2v 


(1 - e C„ 


(18) 


(19) 


Thus 


s 2n 


JX 

[ »,(1 - . V 

r 

s': 


dE 


d s a 2 r 2n (l - e 2 ) 2 (l - e C ) 

(T7T77 “• c * ] “ 


4jt r •'O 


|i a /*2jt 


4jt r 


(1 ■ e C„ 


1, 


(2 - e 2 ) C 2 - 2e C E + 2 e 2 - 1 

d-S ) 2 


dE 


(20) 
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Integration yields 


4 a /»2jt 


4n r 


0 


,3 a 

c e } dE = T^ 

4 r 


(2 + 3 e Z ) 


|ia y* 2jt 


4jt r * / 0 


f * [( 2 - e 2 ) C 2 - 2e C E + 2e 2 - l](l - • C £ ) 


5 a 
4 


2 2 


( 21 ) 


so that the Sun's disturbing function with the short-period terms removed is 

* [¥ » 2 ■«„ * « •“ * * < 2 * 3 * 2 >(! e - ■)] <22) 


4 a 


4 r 


or, substituting for d, $, and C: 
2 


R = 


s 

f c L S L S i c n-A “ C L c i c n-A S a-A ] s 2to 

L S S S s s sj 

PA * K b -\ • v-)1 • ■]} 


+ 15 e 


+ (2 + 3 e 


(23) 


The next step is to remove the medium period effects caused by the motion 
of the Sun. One assumes that the central angle of the Sun 0 g (see Fig. 2) 
varies considerably faster than either co or Cl , and hence R g may be 
averaged with respect to the motion of the Sun. The variations of CO and 
Ci during one revolution of the satellite are given by 
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where the mean orbital rate of the satellite is defined as n = /a 3 and 

the semi-latus rectum of the orbit p is written 

2 

p = a(l - e ) = r p (l + e) (26) 

Defining the mean motion of the Sun n to be 

s 


(27) 

with a g the serai -major axis of the Martian orbit (equivalently, the apparent 

Sun's orbit around Mars), the condition that db and ft are (assumed) much 

smaller than n leads to 
s 



T „ 2 !/ 2 
J 9 

2 e m 

7/2 , 2,2 << "s 

a (1 - e ) 


(28) 


This condition then is satisfied throughout most of the "region of interest." 

Introducing the solar coordinate 0 g and the inclination I g (see 
Fig . 2 ) through 
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c e = c a c l 


s l = s i s e 


C L S A = C I s e 


the function R takes on the somewhat lengthy appearance 
s 

= - ^3 (-%) [ 2 + 3e2) { 1 - l( c a c 20 s + c n + 2 Wi> s 


4 a \r 


2 2 

" s n c i c 26 + s 


2 c V- ffcVc +C V- 

n I J 4 [\i a 2e s in 

\ / 2 2 2 
+ c i s i s o s i ) s 2 e " \ c i c n c i 

3/ 8 v 


2,c i s n c n c i 


+ 2c i s ! c n s i +s i s i V 


+ < c i c n c i +2c i s i c n s i c i +s i s i 


15 2 
4 ® 


2 2 2 2 
c i s n c 2 e +c i s n 


'26 


+ ( c i c n c i +s i s i ) ” c n c 2 e - c n +2c n s n c i> 


(29) 


' 2 ( c i s n c n c i g + c i s i s n s i s ^ s 2 e s _ ( c i c n c i g + s i s i | 

2 2 2 2 
- s n Cl >/ s n c i 8 JS» 

i s n - 2 ( c i s n c n c i s + Wn s i g ) s 2 e s ‘ ( c i c n c i g + s i s i g ) 


2 2 2 _ 
c i s n c 2 e +c ^ 


20 


*(#, ‘Vi ) 


2 2 


2 2 


(30) 


As before, the following trigonometric identities are introduced: 
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C 26 “ C 2 V C 2m 

S 2v S 2o) 


s s s 

s s 


S 2e = S 2v 

’ S 2oj C 2v 

(31) 

s s s 

s s 


and recalling relations 7, applied to the 

Sun: 




A further simplification is introduced here by assuming that the Sun’s 
orbit is circular so that 



The slowly varying disturbing function due to the Sun then becomes: 
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The total slowly varying disturbing function due to J ^ and the Sun per- 
turbation is given thus in Keplerian elements of the satellite and the Sun by 
Eqs . (8) and (34). 


III. LONG PERIOD VARIATIONS IN ECCENTRICITY AND INCLINATION 


In the previous section the slowly varying disturbing function, periodic 
in only Cl and oo, was determined. The long-period changes in the orbit 
elements corresponding to that disturbing function can be readily determined 
by invoking the techniques of Hamiltonian mechanics and canonical transforma- 
tions/ 3 ^ If one introduces the slowly varying Hamiltonian through 




(35) 


where L = (\x a . 

m ' 


then the canonical equations of motion in terras of (slowly varying) Delaunay 
elements are 


L = 


' 51 




Sh 

= 51 




(36) 


where 


iJ = M , g = cd , h = a are 



generalized coordinates 
H = G cos I are generalized momenta 


The set of Eqs. (36) may be integrated if one can find a suitable canonical 
transformation, determined by a function A, such that the transformed Hamil- 
tonian is a function of the new momenta only. Such a transformation, from the 
old variables (L, G, H, i , g, h) to a new set (L', G\ H\ l ' , g' f h’), 
will be given by the function & which is assumed to be expandable in powers 
of the parameter (^ is order Jg): 
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£ = + gG' + hH r + ^(l*, G', H\ i t g t h ) + 


(37) 


Since the canonical form of the equations of motion is preserved and the new 
Hamiltonian is a function of f, G\H* only, it follows that L- c ' H- 
are constants and h* are linear functions of time. The ’old’ 

variables are related to the new by the formulas 


bs 

= bl • 

r d,? 

H ^ 

" = §h 



bs 

= §i7 ■ 

g- = ! L 
Sg 7 

II 

sS 

bi 

Sir 

(38) 


In order to apply the theory to the present problem, the disturbing function 

bS Wrltt6n ln termS ° f «e^unay variables. It turns out to be con- 
venient to write it as a sum of three parts, which follow from Eqs. ( 8 ), (34 ) 

and the definitions (36), and to introduce the perturbing Hamiltonian H • 

P 


- K = R = R +r+r 
P m+s ss m + K sp 

where the terras independent of g and h are 

2 4 


(39) 


R 




4 2 

r J 2 R e / H 2 \ 

and those periodic in g an d h 


(40) 


"s L " 

Sp = 77* K C 2h + a 2 C h + «3 c 2g + «4 C 2h +2 g + « 5 C 2h-2g ' Vh + 2g + « 7 C h-2g] 


( 41 ) 
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The coefficients a ± are given by 


3 2 
a i " 8 S I 



X -»! 
1 2 
G 


■ 8 -!. (-?>*«’ 


i / 2 


15 „2 
a 5 = 16 I 


“* ' 8 V 1 . (‘ " ‘ ?) 


, 1/2 


(42) 


The theory of this study assumes that the function R* dominates R^. If 
the following approximations are used for R m and R gs 

% J 2 R e 

R = 


• = a 3 (l - e 2 ) 3 ^ 


R £ n 2 a 2 (2 + 3e ) 


(43) 
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then 


a « 



111 ^ G 




(2 + 3e ) 


1/5 


(44) 


This inequality holds throughout most of the "region of interest" (see Fig. 1). 

The Hamiltonians associated with the three principal parts of the func- 
tion R are 

m+s 


H - - R 
sp sp 

K - - R 


K = - R 
m m 


K 

s 



- R 

sp 


(45) 


Note that since l is already ignorable we may assume A to be indepen- 
dent of £ so that L = L* = constant. Utilizing Eqs . (38), the perturbing 
Hamiltonian in terms of mixed variables (l 1 , G 1 , H 1 , g, h) becomes 


- - ( dj bS\ _ / SJ bA\ 

*P = H « l L '- G ’ + ST «’ + at) + «ss V L '' G ' + "5F’ H ' + st) 


- ( ^1 ^1 \ 
H sp( L '- G ' + ^T- «' + ^F- h ) 


(46) 


Expanding in Taylor’s series about G' and H f and retaining only one term 

in K and H yields 
ss sp 

H = K (L\ G', H') + H (L 1 , G', H‘) + — -vT + — 
p m 38 ^ <*• ^ 


M (L', G', H’, g, h) 
sp 


(47) 
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In order that K be a function of only (L*, G* , H') one requires that 

(48) 


K (L',G , f H , ,g t h) + — - -g— + — ~ -gj- = 0 
S P ^G' dg ^H‘ dh 


This leads to the condition that 


4 2 


^ R e~ J 2 / ,H' 2 \ S< 1 3 % R e J 2 

-* 3 - 5 ^ 


2L ' G’ 


„ 2 L' 4 


2~ [a L C 2h + °2 C h + a 3 C 2g + a 4 C 2h+2g + a 5 C 2h-2g 


4 u 


" °6 C h+2g + a 7 C h-2g ^ ° 


(49) 


Therefore , choose ^ so that 

3 ,4 n 2 L' 4 

\ = - ~~~4 ' ~ 2 ’ ~ 2~ ^ll S 2h + a 22 S h + a 33 S 2g 

3 J 2 R e 4 


+ a 44 S 2h+2g + a 55 S 2h-2g + a 66 S h+2g + a 77 S h-2g 1 


where the quantities a are to be determined. Substituting into 

Eq. (49) and taking the partial derivatives, the following condition is 
obtained : 

a i C 2h + a 2 C h + °3 C 2g + a 4 C 2h + 2g + 0 5 C 2h-2g " a 6 C h + 2g 

Vt* * (‘ ' 5 jS) <** 


(50) 


+ a_ 


: 33 C 2g + 2a 44 C 2h+2g 2a 55 C 2h-2g 
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+ 2a 66 C h-2g ' 2 ° ! 77 C h-2g^ + 2 G' ^llSh + ®22 C h 

+ 2a 44 C 2h + 2g + 2 «55 C 2h-2g + «66 C h + 2g + «77 C h -2g ) = 0 <») 

Then equating coefficients 



With the Of^'s substituted into (50) the function 4 takes on the 
appearance 
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The coefficients of the trigonometric terms in £ contain six critical 
divisors. An examination of these divisors reveals eleven critical inclina- 
tions which are summarized in Table 1. 

From Eq. (53) the long-period variations in the elements & , g, h, L, G, H 
can be found by appropriate partial differentiations. In this investigation, 
the variations in e and I are of primary importance, and are obtained as 
follows . 

Form the differential of G from its definition (Eq. (36)) 

&G = &[|-i a(l - e 2 )] 1 / 2 
m 


or 


&e 


&G 


\x ae 
m 


G 

2 

L e 




1 


~5g 


Substituting for the variation of 


is given as 


(54) 
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6e = 


2 T . 5 ^» 5 

n L' G* 
s 

6 2 
3 p R e 
m 2 e 


(‘S-M 


"2g+2h 


(•.‘S-s-r-'w-f-y 


h+2g 


/h- 2 h- \ h - 2 « 

( 5 ^ + --\) 


which becomes 


2 2 


5e 


B? (1 + Cj,) 2 

S 

4 (s C 2 , - 2 C r - l) 


(“KH 


( x - 5 c i-) 


2cjd* 


2a’+203’ 


S I t 1 - 0 !.) 2 

s 

4(5 C 2 , + 2 C ir - 1 ) 


(55) 


" 20 ’-ao 1 


s i c i s i- (l - c r> 


S i s c i s s r (l + c r> 

( 5 Cj, - c r - 1) a ' +2u> ' 5 c 2 ' + c It - 1 n ' 


-2co 


(56) 


In order to simplify the equation and the checking of dimensions the parame- 
ters n (mean motion of satellite) and p (defined previously) have been 
introduced. It appears from an inspection of Eq. (56) that at a critical in- 
clination the amplitude of a term in 5e becomes infinite while the frequency 
approaches zero. Actually, as will be shown later, near any critical inclina- 
tion, other than I s= 90°, 5e experiences a finite maximum variation, 
c 

To obtain a corresponding expression for 51 it is convenient to express 
cos I* in terms of H' and G' by using 
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Taking the differential of both sides 


(57) '• 



m v ~ ' I 


To assist the reader in evaluating the equations of variation the partial 
derivatives of ^ with respect to the variables g and h are shown below 
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The variation in I written in canonical variables i 



and in Keplerian elements 
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Sj (1 + Cj,) 2 ^, 


s i (1 - c r ) c r 


+ / \ C 2n'+2CD' / 2 \ 

4 ( 5 c i ’ " c i ' - 7 4 v 5 C I , + 2 C I ' " V 


c 2n , -2o)' 


S I C I S I' C I' (1 + s i C I S I ,C I' (1 " C I ,) 


( 5 c v - c r • x ) 


C f2 , +2o)' + 


( 5 c v + C I- - x ) 


a'-fco' 


2 „2 

s i s i- 


+ T%T c 2a' + s i c x s i ,c n’ 


s 2 (i + c v r 


s i (1 " c i^ 




4(5 Cj, - 2 C 


s x Cj (l + c r )s It 


C 2 n'+ 2 u>' + 


4(5 C 2 , + 2 C,. - l) 
s ! C I - c r )s i- 


C 2U ' -20)' 


+ “7 9 \ C n , +2a i t + / 2 \ 

2 ( 5 c v ■ C I' ' V 2 \ C I’ + C I' “ 1 j 


'n'-au' 


(61) 


Note here that the inclination I c = 90° is critical for SI but does not 
appear in the expression for be. 


IV. LONG-PERIOD BEHAVIOR NEAR CRITICAL INCLINATION 

The technique outlined in the previous section for obtaining long-period 
behavior of the elements proves unsatisfactory when the inclination is near 
one of the critical values (see Table l). In this section a suitable analysis 
will be conducted to determine the behavior of the satellite's motion while in 
the neighborhood of a general critical inclination. From Eq. (41 ) one can 
write the periodic part of the Hamiltonian in the form 
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n 2 L' 4 
s 


sp 


4 


1 


i=-2 


7 ± cos (2g» 


ih» 


^ ^ cos jh' 


(62) 


J-1,2 


The critical inclinations are just those at which one of the trigono- 
metric arguments in (62) has a zero rate due to J 2 . Thus the terms of 
cos jh' yield a critical inclination of 90°; however, since they do not con- 
tain the variable g' they will contribute nothing to the variation of e'. 
Of course, these terms will contribute to the long-period behavior of the 
inclination, but since the primary emphasis in this investigation is centered 
on the radius of pericenter (equivalently e f , since 5r^ = -a'Se'), these 
terms are not of direct interest here. Thus in the following analysis these 
terms will also be considered to have been removed by a suitable generating 
function, so that the total Hamiltonian of interest will be 


^ 7 i cos (2g» + ih' ) (63) 


2 4 2 

n L’ 
s 


m i=-2 


Specifically, the coefficients y^ are given in terms of the of 


Eq. (42) as 


7 o = a 3 


>1 = ^6 
7 2 = a 4 


7 -l = “7 
7 -2 = «5 


(64) 


In order to analyze the motion near the critical inclination defined by 
2g* + ih* = 0, a determining function £ is introduced which is taken to 
be £ with the i term of interest absent. Thus £* will "remove" all the 
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periodic terms from the Hamiltonian except the i term, so that H p ( = " R m+g ) 


now reads 


*2 T ,4 
n L’ 


R = -H = R + R + 7 T- 7. cos (2g ' + ih') 

m+s p ss m . 2 'i 


4 V- 


(65)' 


Introduce now a contragredient transformation to new variables 


Ik" 1 f 1 

i/ 2 ] iv ■ 

1 h” J [o 

1 J 

L-J 

n r x 

°1 

- G t- 

Lh"J L-i/2 

J 

H'_ 


Then the Hamiltonian is of the form 


R* = R^ (G", H m , L') + R* (G’\ H", L’) + 
m+s ss v m 


( 66 ) 


(67) 


n 2 L ' 4 + 

— — g- 7 cos 2 g" ( 68 ) 

4 


Note that L ' and H" are constants and, since the Hamiltonian itself 

is a constant in the absence of explicit time terms it follows that R m+S 

also a constant. Define now G" as that function of L' and H" for which 

c 


■ ° 


and expand R about (G" - G" ): 

m+s c 


R* SR* (G", H", L') + R*(G", H", L') + [^7 R* 1 (G" - G”) 

m+s ss c * m e LoG" ss J G , t=G »t c 


L^g" 


(g m -g m ) 


JG"=G" 
c 


n 2 L ' 4 


cos 2g M (69) 


G =G 
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Ignoring the fluctuations in y by comparison with its value at G" is not 

1 c 
valid for small eccentricity (see Appendix). Using the definitions 



(70) 


equation (69) may be written as 

V G " _ G c )2 + ' G”) + ® 3 + f cos 2g" = constant (7l) 


Completing the square and suitable rearrangement yields the equation 

\ , S 2 f f 

1+ + « cos 2g n = K = constant (72) 

/ 2^ H"J H" 

which is now recognized as being in the form of the simple pendulum equation. 
The solution of this problem is well known and results in a plot of phase space 
contours which describe the motion of the pendulum. In the present problem, 
the variation of G ,, /h" with the angle g" follows from 



91 

H" 




f 

® x H " 2 


-.1/2 


cos 2g" 


2® x H" 


(73) 
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Note that g" can librate about n/2 or 0 according to whether > 0 

or <0. As in the pendulum problem, the separatrix is determined by setting 
K = f/SB^H"^ and varying g". The motion bounded by this contour represents 
a libration of C'/h" about some average (G ,, /H” ) ave » and as can be seen in 
Fig. 3, the maximum fluctuation of this variable follows a contour just 
inside the separatrix and is equal to 2 yfi K. 

Recall now that it was assumed earlier that R gg « R^, so that 
can be taken as 




(74) 


J R 2 / 

3 m 2 e L 

1 2 3 5 V 


cos 1 1 - lOi cos I ' + 




which, written explicitly in terms of the old primed variables, is 

4 _2 

^m_ 

L* ~G ’ 

The maximum variation in G"/h" written in terms of the original orbital 
elements is 

W 2 


(75) 


6 © = 2 = 2 


2f 


S H 


i . 2 


(76) 


from which the maximum variation in the inclination is found to be 


BI 


1 max 


. , , , ,2,3/4 ,5/2 

2 cos 1 1 -i e * (l-e ) ' n a' ' 

1 c 1 s 

5 

R sin T 

2 \ 

/ 2 . i „ \ | 

e c 

V 2 I ' 30 COS I c" 101 COS I c + T + 4 ^I 


|l/2 


(77) 


since 


■©- 


sin 1 1 


51=2 


(“■ ■; 


-I) 


2f 


IB H n " 


1/2 
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* * 
with 7 defined by y 


15e 


,2 'l 


7 t (I 




The theory developed near 1^ is valid if 61 is small enough so that 
I' does not get too near any other neighboring r . To prevent this occur- 
rence, an upper limit is assigned to the maximum variation 1 51 1 such that 

1 1 max 


I Si I 




(78) 


where l£ is the 1 ^ under investigation and is the nearest neigh- 

boring one . ^ 

The corresponding maximum variation in the eccentricity is found from 
the constancy of H" i.e., 5H" = 0 which gives 


6e 1 


2(1 - e* ) sin I* 


i 1 1 2 cos 1 1 - i I 


|8I| 


(79) 


l 5e ’l, 


2n s a 5 / 2 (l-e' 2 ) 3 / 2 


5 7, 


U J l(30 cos 2 X’ - 101 cos I 1 + - 4)1 

ra 2 1 \ c c 2 /I 


Now define the constant by 


K = 2n 
i s 


R 3 / 2 


5 7 . 


pi J I /cos 2 I 1 - iOi cos I 1 + 1 4)1 

m 2 1 \ c c 2 / 1 


1/2 


1/2 


(80) 


(81) 


so that 


|8e’| = K. (l 

1 1 max i 


. 2 ) 7/4 


(t) 


5/2 


(82) 
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Then 


61 


can be written in terms of 


SI 


| 2 cos 1^ - i I 
2 sin I* 


K ± e'(l-e' 2 ) 3/4 


(t)' 


5/2 


(83) 


Invoking then the condition on a maximum allowable | AI l 11 , ax ( E< 1- ( 78 )) 
one can write an upper limit a' as a function of e* consistent with the 
analysis as 


kL 

R 


M 

^ 1 ' max 

f 2 sin i£ A 

■ ij 

n 

y/l 

| 2 cos I* - i | K . 
U c 1 i 

l_e ' (l-e ' 2 ) 3 ^J j 


2/5 


(84) 


or, defining constants by 


2 sin 1 ' A 


I 2 cos I ' - i K 4 


2/5 


then 


(r ) 


P max = B e ’~ 2 / 5 ( 1 +e 1 )" 3 / 10 (l-e 1 ) +7 / 10 
R i 


(85) 


( 86 ) 


The constants B j , and K ^ for the different values of i(-2 < i < 2) 
corresponding to the 10 different critical inclinations are listed in Table 2. 
The bounds (86) on the applicability of this analysis are shown in Fig. 1. 

Finally, knowing the maximum fluctuations in e’, the maximum variation 
in the radius of pericenter is determined as 


l Br , 


p 1 max 


a 1 


r 7 / 2 

_P 

i ^ 



(87) 


This maximum variation of the radius of pericenter with e' near the different 
critical inclinations is shown in Fig. 4 for r^ = 4000 and 6000 km. For r^ 
above the theoretical limits of Fig. 1, the corresponding portion of the curve 
is dotted in Fig. 4 for a specific index i, the theory being open to ques- 
tion since the fluctuations in I 1 may overlap a neighboring critical 
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inclination. In case of overlap a simple analysis of |&e| is not 

1 1 max 

possible . 

This investigation was purposely limited to large values of e' due to 
the mission requirements imposed on it; in fact, the analysis of Section IV 
was invalid for e' close to zero. For e f close to zero, however, the 
maximum variation in e 1 agrees with that found in this analysis; this agree- 
ment is demonstrated in the Appendix. 

APPENDIX 

In this section, the eccentricity is assumed to be of order j and 

it is assumed that the inclination is near one of the critical inclinations 
I£ other than 90° . 

Recall that in Section IV it was found that the Hamiltonian could be 
written in the form 

a i (G" - G^') 2 + S 2 (G" - G") + S 3 + f cos 2g” = constant (A.l) 

Then from the constancy of H” 

H” = H' - | G' = const 


it follows that 


A - ®' 2 ( cos i' - 1) 


const 


in (l-e* ) + 2 in |cos I* - i) = const 

e' 2 = 2 in |cos I* - — j + const + 0(«J 2 ) 

Next , expand about I » = I * : 

c 


(A. 2) 


(A. 3) 
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i' 2 = 2 




A V Sl11 1 ' 

(cos i; - 2 ) + — “T e -t (i '" i c )+ ••• 
cos r - 2 


„ 2 sin I' 2 

e- = S_.(I'-I') +H, + 0(J ) 

cos V-- 


where H is a small constant of order J_ if the initial values of 
1 ^ 

and I'-I* are 0(j _). Therefore 
c * 


...... ti '1 'J l 

1 c 2 sin I' [ 


From the constancy of H" 

&H" = SG" (cos r - sin r 61' =0 


so that 


G" sin I* 

&G" = — j 61 1 

cos v - 2 


Substituting 61* = I * -1^ from (A. 5) 



the Hamiltonian now takes the form 





, 2 



2 G" r 

+ *2 T [ e 


2 


H J 


2 # 

+ e* f cos 2g" = constant 


where 


f/e 


2 


2 

^This form of f* does not, in fact, contain a l/e' factor since 
of definitions 70, 64, and 42 reveals that f has an e' factor. 


(A.4) ' * 
e ' 2 

(A. 5) 


(A.6) 

(A. 7) 
(A. 8) 

inspection 
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Combining powers of e 1 , Eq. (A.8) is rewritten as 


Sj©’ 4 + (® 2 - 2 e' 2 + f*e' 2 cos 2g" = constant (A. 9) 

G " 2 _ G" 

where 3, = S, , 3, = 3 g ^ . 

The interplay between g" and e f (and hence I') may be more easily 
visualized if one introduces the coordinates 

| - e f cos g" 

T] = e' sin g" (A. 10) 

so that Eq. (A.9) becomes 

+ n ) + (® 2 - 2 S 1 H 1 )(| 2 + T) 2 ) + f*(| 2 - T) 2 ) =x (A . 11 ) 

The equilibrium points of the contours in the |-T| plane are found from 

|f = s [ 2 + T! 2 ) + (s 2 - 2 + f*J = 0 (A. 12) 

= nj 2 I^! 2 + i, 2 ) + (» 2 - 2 SjH^ + f*J = 0 (A. 13) 

Solution of these equations yields equilibrium points at the origin of the 
I > ^ plane and on the £ and rj axes at 


and 


respectively. 


( 0 ,t! 2 ) = (o, 


(S 2 ,o) = o, 


(s 2 + f*) 

2S. 


|l/2\ 


* — 
f - ®o 


H 1 + 


2 3, 


1/2 


- 0 
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From (A.ll) one can see that the contours are symmetric about both axes; de- 
pending on the value of a variety of contours are possible about the 

equilibrium points , c.f., Fig. 5. 

The "separatrices" between different types of contours are shown by 
dotted lines. As increases, (a) “► (c) (e) or (b) (d) -► (f) and 

the separatrices, after appearing in (c) or (d), grow and then change to cir- 
cular form in (e) or (f). In (g) and (h) are shown their transition forms 
corresponding say, to an< * ® s ^ own 

slightly larger. 

On contours outside the circle-pair, and for f < 0 for example, 

from (A.ll) one finds that e (corresponding to y = 0) is given by 

max 

I 2 = J P~+ y/p+k (A. 14) 

max 


where 


x/7 = 


*2 + f 


28, 


while e (corresponding to £ = 0) is 

min 

e = J~Q + \/Q+k ( A . 15 ) 


Then 

(8e) (e + e 4 ) = e 2 - e 2 . = Jv-yfc + - 7 == "4= ( A - 16 ) 

v max max min max min v v y p+k + 

which is a decreasing function of k while e + e is an increasing 

° max min 

function of k. It follows that (Be) decreases as k increases, and 

max 

hence is smaller outside the circle-pair than on it. 

Thus it is apparent in Fig. 6 that the largest fluctuation, Be, in eccen 
tricity occurs on a contour just inside a separatrix circle. It will be shown, 


- f 


where ^/q = H — 

258, 
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moreover, that this ( &e ) fflax remains constant as 1^ increases further in 
spite of the growth of the radius of the separatrix circles. The overall 
maximum 5e is thus just the diameter of the circles in (g) or (h). To 
evaluate this, note that for sufficiently large H, the family of contours 
(A.ll) includes the circle-pair 

[(8 - y ) 2 * n 2 - p 2 ][(i + 7 ) 2 + n 2 

or 

[5 2 + (rj - 7 ) 2 - P 2 ][| 2 + (n + 7 ) 2 



But ( &€ ) max “ (p + y) - (p “ 7 ) = 2 7» independent of p. Thus 


P 2 ] = 0, < 0 (A. 17) 

2 f* 

P ] = 0, — > 0 (A. 18) 

IB. 


(Se 


= 2 




(A. 19) 


in agreement (for small e') with Eqs. (77) and (79) of the large e 1 theory, 
which yield 


(&e) 



c 



(A. 20) 
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TABLE I 

CRITICAL INCLINATIONS 










HIGH ECCENTRICITY ORBITS ABOUT MARS 
REFERENCES 

1. Garfinkel, B., Astronomical Journal , 1960, Vol . 65. 

2. Tlsserand, F., Traite de Mecanique Celeste , 1889, Vol. 1 

3. Brouwer, D. , Astronomical Journal , 1959, Vol. 64. 

National Aeronautics and Space Administration 
Electronics Research Center 
Cambridge, Massachusetts, December 1966 
129-04-04-07 


209 


TRAJECTORY ANALYSIS AND GUIDANCE THEORY 



0.4 0.5 0.6 0.7 0.8 0.9 

e 


BOUNDARY CURVES FOR VALIDITY OF CRITICAL INCLINATION ANALYSIS 

Figure 1 
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CONTOURS FOR H, SLIGHTLY GREATER THAN H, c 
Figure 6 
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